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Recently  a  number  of  new  techniques  for  constructing 
nonlinear,  second-order  accurate,  high-resolution  shock- 
capturing  schemes  for  systems  of  hyperbolic  conservation 
laws  have  been  developed  in  order  to  simulate  complex  flow 
fields  more  accurately  and  efficiently.   These  schemes, 
called  total  variation  diminishing  (TVD) ,  are  free  from 
generating  spurious  oscillations  across  the  shocks  and 
contact  discontinuities,  and  can  converge  to  a  physically 
realizable  solution  with  the  aid  of  entropy  correction. 

The  extension  work  has  been  conducted  to  extend  the 
TVD  schemes  developed  for  two-dimensional  planar  flows  to 


three-dimensional  flows.   On  the  application  side,  a  thin- 
layer  Navier-Stokes  flow  code  based  on  the  conventional 
Beam-Warming  scheme  has  been  modified  to  incorporate  the 
TVD  schemes  in  the  solution  process.   A  variable  time  step 
procedure  has  also  been  implemented  into  the  code  for  more 
efficient  computation  of  the  steady  state  flows.   The 
transonic  turbulent  flow  over  an  axisymmetric  projectile  is 
adopted  as  the  test  problem  to  assess  the  performance  of 
three  variants  of  the  TVD  scheme,  two  upwind  schemes  with 
flux  limiters  proposed  by  Van  Leer  and  Roe,  and  one 
symmetric  scheme.   The  results  obtained  show  that  all  three 
variants  can  accurately  resolve  the  shock  and  flow  profiles 
with  fewer  grid  points  than  the  Beam-Warming  scheme. 
Furthermore,  although  the  computing  cost  of  the  TVD  schemes 
is  higher  than  the  Beam-Warming  scheme  for  each  time  step, 
the  higher  convergence  rate  of  these  schemes  makes  them 
computationally  more  efficient  overall.   The  combinations 
of  the  high  accuracy,  good  robustness,  and  improved 
computational  efficiency  offered  by  the  TVD  schemes  make 
them  extremely  attractive  for  computing  flow  with  shocks. 
Among  the  three  variants,  it  is  found  that  the  convergence 
rate  of  the  flux  limiter  proposed  by  Van  Leer  for  the 
upwind  TVD  scheme  is  least  affected  by  the  incoming  flow 
Mach  numbers. 
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CHAPTER  I 
INTRODUCTION 

Computational  fluid  dynamics  (CFD)  has  recently  become 
an  increasingly  powerful  tool  in  the  aerodynamic  design  of 
aerospace  systems.   A  number  of  CFD  codes  have  been 
developed  to  simulate  the  high-speed  flow  problem.   Routine 
calculations  can  now  be  performed  to  partially  substitute 
for  the  wind  tunnel  measurements.   Consequently,  designs 
can  now  be  analyzed  and  improved  in  shorter  time  and  with 
lower  cost.   However,  despite  these  fast  advancements,  much 
remains  to  be  done  to  improve  both  the  accuracy  and 
efficiency  of  the  numerical  flow  calculations.   One  major 
topic  particularly  relevant  to  aerospace  applications  is 
the  treatment  of  shocks. 

In  general  there  are  three  major  difficulties  in 
computing  flow  with  shocks.   These  difficulties  are  listed 
as  follows: 

(i)  Schemes  that  are  second  (or  higher)  order 
accurate,  such  as  the  well-known  methods  developed  by  Lax 
and  Wendroff  [1],  MacCormack  [2],  and  Beam  and  Warming  [3], 
may  produce  spurious  oscillations  across  the 
discontinuities.   First  order  schemes,  on  the  other  hand, 
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suffer  from  severe  numerical  diffusion.   They  cannot 
provide  accurate  solutions  for  complicated  flow  fields 
unless  an  excessively  fine  grid  distribution  is  adopted. 

(ii)  The  numerical  schemes  may  develop  nonlinear 
instabilities  when  discontinuities  are  encountered. 

(iii)  The  schemes  may  yield  nonphysical  solutions. 
It  is  well  known  that  these  classical  schemes  are  not 
robust  and  accurate  enough  for  flow  containing 
discontinuities . 

In  order  to  suppress  high-freguency  numerical  errors 
of  second  or  high  order  schemes  associated  with  the  shock 
calculation,  the  most  common  procedure  is  to  add  artificial 
dissipation  terms.   For  example,  in  the  widely  used  Beam- 
Warming  scheme,  a  fourth-order  explicit  dissipation  is  used 
to  control  non-linear  instabilities,  whereas  the  second- 
order  implicit  dissipation  is  built-in  to  relax  the 
stability  bound  from  the  explicit  fourth-order  dissipation. 
All  these  numerical  dissipation  terms  contain  adjustable 
parameters.   For  inviscid  steady-state  airfoil  calculation, 
the  use  of  the  above  type  of  dissipation  terms  may  still 
produce  spurious  oscillations  near  shocks  [4].   Such 
oscillations  are  not  only  undesirable  but  often  lead  to  a 
breakdown  in  stability  of  the  numerical  scheme.   This  is 
the  major  drawback  of  the  Beam-Warming  scheme  and  other 
schemes . 
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In  order  to  overcome  these  difficulties,  a  number  of 
new  techniques  for  construction  of  nonlinear,  second-order 
accurate,  high-resolution  schemes  for  hyperbolic 
conservation  laws  have  been  developed  in  recent  years  in 
order  to  simulate  complex  flow  fields  more  accurately  and 
efficiently  [5].   All  these  schemes  share  the  following 
attractive  properties: 

(1)  The  schemes  do  not  generate  spurious  oscillations. 
They  are  total  variation  diminishing  in  the  nonlinear 
scalar  case  and  the  constant  coefficient  system  case  [5]. 

(2)  They  are  consistent  with  the  conservation  law  form 
and  entropy  inequality  [6].   This  property  guarantees  that 
the  numerical  weak  solutions  of  the  hyperbolic  conservation 
law  are  physically  realizable. 

In  the  following,  a  brief  review  is  conducted  to  trace 
the  original  ideas  that  the  modern  TVD  schemes  are  based 
on,  and  to  present  several  variants  of  the  TVD  scheme 
proposed  in  the  literature  which  have  proven  to  be  useful. 

CIR  Method 

The  earliest  method  for  gas-dynamic  equations  in 
characteristic  form  was  proposed  by  Courant  et  al.  [7]. 
Their  procedure,  sometimes  called  the  CIR  (Courant- 
Issacson-Rees)  method  is  to  trace  back  from  (iAx,  (n+l)At) 
all  three  characteristic  paths.   Since  the  problem  is 
nonlinear,  the  directions  of  these  paths  are  not  known 


exactly,  but  to  a  first  approximation  they  can  be  taken 
equal  to  their  known  directions  at  (iAx,  nAt) .   Then  each 
characteristic  equation  is  solved  by  using  interpolated 
data  at  time  nAt,  in  the  interval  to  the  left  of  i  for 
characteristics  with  positive  speed,  and  in  the  interval  to 
the  right  of  i  for  characteristics  with  negative  speed. 
This  method  has  two  principal  drawback:  because  it  does  not 
deal  with  the  integral  form  (conservative  form)  of  the 
equations,  it  cannot  convey  a  shock  wave  with  the  proper 
speed;  also,  the  resulting  method  is  only  first-order 
accurate. 

Godunov's  scheme 

In  1959,  Godunov  [8]  proposed  the  idea  of  advancing 
the  solution  to  the  next  time-level  by  solving  a  set  of 
Riemann  problems.  It  is  noted  that  the  Riemann  problem  for 
any  hyperbolic  system  of  conservation  laws  arises  if 
initial  data  are  prescribed  as  two  semi-infinite  states  (u 
=  uL  for  x  <  0,  u=  uR  for  x  >  0  )   For  one-dimensional 
Euler  equations,  the  solutions  consist  of  three  waves  (u, 
u±c)  which  represent  the  propagation  speed  of  the  contact 
discontinuity,  expansion  wave,  and  shock  respectively. 

To  describe  Godunov's  use  of  the  Riemann  problem,  let 
uj+l/2  °e  the  average  state  over  (j±l/2)Ax  at  (nAt).   He 
then  replaced  the  data  by  an  approximate  distribution  in 
which  the  state  inside  each  interval  is  uniform  and  equal 


5 
to  u j .   For  each  interface  (j+l/2)Ax  one  can  solve  the 
Riemann  problem  with  ur  =  uj+^  and  ul  =  u j .   This  gives  an 
exact  solution  Uj+1/2  to  tne  approximate  problem,  assuming 
that  At  is  small  enough  that  the  waves  do  not  interact  from 
neighboring  interfaces.   The  solution  at  (n+l)At  can  again 
be  approximated  by  a  piecewise  uniform  distribution,  and 
then  the  process  can  be  repeated. 

The  advantage  of  Godunov's  scheme  is  the  clear 
physical  picture  of  interaction  involved,  and  the  absence 
of  any  arbitrary  parameter.   Even  though  Godunov's  scheme 
is  monotonicity-preserving  and  produces  good  steady  shock 
profile,  the  main  drawback  of  Godunov's  scheme  are  that  it 
is  cumbersome  to  compute  and  that  is  involves  the  high  cost 
to  solve  the  Riemann  problem  exactly.   Also  this  scheme  is 
only  first  -order  accurate. 

Flux  Correct  Transport  (FCT)  Scheme 

The  flux  Correct  Transport  (FCT)  scheme  was  first 
proposed  by  Boris  and  Books  in  1973  [9],  and  later  extended 
and  refined  by  Zalesak  [10].   The  FCT  scheme  may  be 
considered  as  a  two-step  method,  a  lower  order  transport 
stage  corrected  by  a  higher  order  antidif fusive  flux  in 
such  a  way  as  to  prevent  the  occurrence  of  spurious 
oscillations.   In  general,  The  FCT  scheme  performs  well 
compared  to  other  schemes.   However,  as  Zalesak  has  pointed 
out,  peaked  data  clipping  can  occur,  and  he  proposed  a 
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modified  flux  limiter  which  eliminates  the  clipping  and  is 
easy  to  extend  to  multidimensions. 

Harten's  Scheme 

In  1983,  Harten  [5]  developed  a  simple  algebraic 
criterion  for  constructing  a  first-order  accurate  TVD 
scheme.  He  then  added  a  limited  anti-diffusive  flux  to  a 
first  order  accurate  entropy-satisfying  scheme  to  obtain 
the  second-order  accurate  TVD  scheme.   This  technique  is 
referred  as  the  'modified  flux  approach'.   The  modified 
flux  is  chosen  so  that  the  scheme  is  second-order  at  a 
smooth  region  and  will  switch  itself  into  first-order 
accuracy  at  the  vicinity  of  points  of  extrema.   In  [11]  he 
also  has  extended  the  class  of  explicit  TVD  schemes  to  a 
more  general  category  which  includes  a  one-parameter  family 
of  implicit  second-order  TVD  schemes.   More  recently 
Harten's  TVD  scheme  was  modified  and  generalized  by  Yee 
[12]  and  implemented  to  solve  the  two-dimensional  Euler 
equations  of  gas  dynamics  for  the  airfoil  problem.   The 
numerical  results  indicate  that  Harten's  scheme  is  both 
robust  and  accurate. 

MUSCL  (   Monotonic  Upwind  Schemes  of  Conservation  Laws  ) 
Schemes 

B.  van  Leer  [13]  adopted  an  approach  developed  by 
Godunov  (1959),  who  used  a  control  volume  formation  to 
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ensure  that  the  scheme  would  be  conservative  and  therefore 
dealt  with  averaged  values  instead  of  nodal  values,  Van 
Leer  also  observed  that  it  is  beneficial  to  work  with 
averaged  gradients  in  addition  to  averaged  values,  i.e., 
one  can  obtain  second-order  accuracy  in  Godunov's  scheme  by 
replacing  the  piecewise-constant  initial  data  of  the 
Riemann  Problem  with  piecewise-linear  initial  data.   The 
slope  of  the  piecewise-linear  initial  data  is  selected  so 
that  spurious  oscillation  will  be  prevented.   His  result 
indicated  that  one  does  not  have  to  resort  to  a  Godunov  - 
type  of  Riemann  solver  in  order  to  obtain  good  shock 
resolution.   Colella  and  Woodward  [14]  introduced  the 
piecewise  parabolic  method  (PPM) ;  They  further  refined  van 
Leer's  idea  by  using  piecewise  parabolic  initial  data,  and 
the  results  they  obtained  are  remarkably  good. 


Roe's  Scheme 

Roe  [15,16,17]  has  developed  a  numerical  fluctuation 
approach  toward  schemes  for  the  numerical  solution  of  both 
scalar  hyperbolic  conservation  law  and  hyperbolic  systems 
of  conservation  laws.   His  approach  combines  a  novel 
interpretation  of  difference  schemes  with  higher  order 
accuracy  and  monotonicity  preservation  and  also 
incorporates  these  ideas  into  an  approximate  Riemann  solver 
for  hyperbolic  systems.  Roe's  approximate  Riemann  solver 
requires  that  there  exists  an  average  function  A(uR,uL) 
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which  is  a  nonlinear  function  of  up  and  ul.   The  average 
function  is  constructed  by  the  use  of  flux  limiters.   Using 
this  approach,  Roe  has  obtained  solutions  with  a  crisp 
shock,  free  from  the  spurious  oscillations  inherent  to 
classical  shock-capturing  schemes.   Unfortunately,  Roe's 
scheme  may  admit  nonphysical  (entropy  violating)  solutions, 
i.e.,  expansion  shocks.   Various  remedies  have  been 
proposed  [18],  [19]/  [20]. 

Osher's  Scheme 

Osher's  scheme  [21]  is  based  on  an  approximate  Riemann 
solver,  but  compression  and  rarefaction  waves  are  used  to 
approximate  shocks.   This  leads  to  a  simpler  and  smoother 
algorithm.   The  numerical  flux  functions,  which  are  at 
least  continuously  differentiable,  are  written  in  closed 
form  and  include  various  switches  which  make  them  upwind. 
Also  the  limit  solutions  satisfy  the  entropy  condition.   It 
is  observed  that  Osher's  scheme  requires  more  operations 
than  Harten's  or  van  Leer's  [22]. 

Symmetric  TVD  Scheme 

In  1984,  Davis  [23]  showed  that  a  certain  class  of  TVD 
schemes  can  be  interpreted  as  a  Lax-Wendroff  scheme  plus 
and  upwind  weighted  conservative  numerical  dissipation 
term.   He  then  simplified  the  scheme  by  eliminating  the 
upwind  weight  of  this  numerical  dissipation  term  and  also 
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ensured  that  the  simplified  scheme  still  possessed  the  TVD 
properties.   Shortly  after  that,  Roe  [24]  reformulated 
Davis's  scheme  in  such  a  way  that  the  resulting  scheme  was 
easier  to  analyze  and  included  a  class  of  TVD  schemes  not 
obtained  by  Davis.   More  recently  Yee  [25]  generalized 
Roe • s  schemes  to  a  one-parameter  family  of  second-order 
explicit  and  implicit  TVD  schemes.   Therefore,  the 
formulation  of  Roe  and  Davis ' s  scheme  can  simply  be  viewed 
as  special  cases  of  Yee's  explicit  symmetric  TVD  schemes. 
Also  the  main  advantages  of  Yee's  formulation  are  that 
stiff  problems  can  be  solved  more  efficiently  by  using 
implicit  method  and  that  steady-state  solutions  are 
independent  of  the  time-step. 

Present  Contributions 

In  the  present  study,  the  viewpoint  of  interpreting 
the  TVD  scheme  as  a  modified  flux  is  adopted  in  order  to 
make  the  task  of  practical  implementation  more  tractable. 
One  can  modify  a  standard  three-point  central  difference 
code  (  e.g.  Beam-Warming  algorithm  )  by  simply  changing  the 
conventional  dissipation  terms  into  the  one  designed  for 
TVD  schemes.   Hence  the  only  difference  in  computation  is 
that  the  TVD  scheme  reguires  more  elaborate  dissipation 
terms  or  flux  limiters  for  the  explicit  and  implicit 
operators.   In  order  to  gain  more  insight  into  the  TVD 
scheme  for  the  complex  flow  problems,  one  objective  of  this 
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dissertation  is  to  examine  how  these  various  flux  limiters 
affect  the  accuracy  as  well  as  efficiency  of  the  numerical 
solutions. 

An  extended  derivation  of  the  TVD  schemes  in  three- 
dimensional  generalized  curvilinear  coordinates  will  be 
presented.   The  formula  will  then  be  applied  to  actual  flow 
computations.   A  time-dependent  axisymmetric  Navier-Stokes 
code  developed  at  the  U.S.  Army  Ballistic  Research 
Laboratory  [26]  has  been  modified  to  utilize  more 
sophisticated  numerical  dissipation  models  based  on  the  TVD 
schemes.   Furthermore,  a  spatially  varying  time-step 
procedure  [4]  is  also  incorporated  into  this  new  version  of 
the  code  to  speed-up  the  convergence  rate  for  steady-state 
flow  computations.   The  transonic  turbulent  flow  over  the 
axisymmetric  secant-ogive-cylinder-boattailed  (SOCBT) 
projectile  with  sting  at  zero  angle  of  attack  is  used  as 
the  test  problem  to  assess  the  performance  of  the  TVD 
schemes  in  terms  of  numerical  accuracy,  robustness  and 
convergence  rate. 

The  transonic  turbulent  flow  over  an  axisymmetric 
projectile  problems  at  three  different  Mach  numbers  has 
been  computed  on  many  different  grid  sizes  with  both  the 
Beam-Warming  scheme  and  three  variants  of  the  TVD  scheme. 
A  systematic  assessment  will  be  presented  to  appraise  the 
relative  strength  and  weakness  among  these  schemes. 
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Numerical  results  will  also  be  compared  with  the 
experimental  data  obtained  by  Kayser  et  al  [27]. 

The  remaining  dissertation  is  organized  as  follows. 
The  basic  properties  of  hyperbolic  conservation  laws  and 
TVD  schemes  are  first  reviewed  in  Chapter  II.   Then  the 
detailed  formulation  of  upwind  and  symmetric  TVD  schemes 
for  nonlinear  scalar  hyperbolic  eguations  and  system  of 
hyperbolic  eguations  with  constant  coefficients  are 
described  in  Chapter  III.   In  Chapter  IV,  A  brief  survey  of 
the  mechanism  of  flux  limiters,  which  play  the  central  role 
in  all  TVD  schemes,  will  be  presented.   In  Chapter  V,  The 
derivation  is  presented  which  extends  the  TVD  schemes  in 
two-dimensions  to  three-dimensions  with  generalized 
curvilinear  coordinates.   In  Chapter  VI,  a  thin-layer 
Navier-Stokes  solver  along  with  some  descriptions  of  the 
related  physical  modeling  and  numerical  procedure  aspects 
are  presented.   Results  of  numerical  computation  and 
detailed  discussions  are  given  in  Chapter  VII,  which  is 
followed  by  a  summary  and  concluding  remarks  in  Chapter 
VIII. 


CHAPTER  II 
PRELIMINARIES 
In  this  chapter  we  give  a  summary  of  important  notions  in 
the  theory  of  hyperbolic  conservation  laws  and  conservative 
finite  difference  methods.  We  also  briefly  review  the  theory 
of  total  variation-diminishing  (TVD)  finite  difference 
schemes . 

2.1  Hyperbolic  Conservation  Laws 

In  this  section,  we  will  consider  numerical 
approximations  of  the  initial  value  problem  for  nonlinear 
hyperbolic  systems  of  conservation  laws: 

(2.1)  Ut  +  [f(u)]x  =0  -co  <  X  <  co, 

u(x,0)  =  u0(x) , 

where  u(x,t)  is  a  column  vector  of  m  unknowns,  the  flux 
function  f(u)  is  a  vector  of  m  components,  and  u0(x)  is  the 
initial  data.   Equation  (2.1)  can  also  be  written  in  a 
matrix  form  as: 

(2.2)  ut  +  A(u)ux  =  0,    A(u)  =  af/au. 
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Equation  (2.1)  is  called  hyperbolic  if  the  Jacobian  matrix 
A(u)  has  real  and  distinct  eigenvalues 


Kl   <  K2   <  <  * 


m 


Here  k^,    i=l, — ,m  ,  are  eigenvalues  of  the  matrix  A(u) . 
It  is  well  known  that  the  solution  of  Eq. (2.1)  may 
develop  discontinuities  even  when  the  initial  data  are 
smooth  [11].   Because  of  this,  we  seek  a  weak  solution 
which  satisfies  Eq.(2.1)  in  an  integral  sense,  i.e., 


(2.3) 


0      J 


[u  <j>t  -  </>x   f  (u)  ]dx  dt  + 


[u0(x)^]dx  =  0, 


where  <f>(x,t)    is  a  continuously  differentiable  test 
function. 

An  immediate  consequence  of  the  integral  relation 
Eq.(2.3)  is  that  every  piecewise  continuous  weak  solution 
must  satisfy  a  jump  condition,  which  is  known  as  the 
Rankine-Hugoniot  condition  in  gas  dynamics,  across  the 
discontinuity,  i.e., 

(2.4)      s  [uj.]  =  [fi]        i=  i,...,m 

Here  the  brackets  [  ]  denote  the  jump  across  the 
discontinuity,  and  s  denotes  the  speed  of  the 
discontinuity . 
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It  is  also  well  known  that  weak  solutions  of 
hyperbolic  conservation  laws  are  not  uniquely  determined  by 
their  initial  values.   A  simple  example  of  such  a 
discontinuous  solution  of  Eq.(2.1)  for  a  scalar 
conservation  law  with  f (u)  =  1/2  u2  is 


ui(x,t)  =  - 


The  discontinuity  is  across  the  line  x  =  t/2,  which 
propagates  with  speed  s  =  1/2;  across  the  discontinuity  [u] 
=  1,  [f]  =1/2,  so  the  jump  condition  of  Eq. (2.4)  is 
satisfied. 

The  function 


1 

for 

X 

< 

t/2 

0 

for 

X 

> 

t/2 

u2( 


x.t,  .  {  ; 


for  x  <  t/2 
for  x  >  t/2 


again  is  discontinuous  solution  across  the  line  x  =  t/2, 
and  [u]  =  -1,  [f]  =  -1/2,  so  that  the  relation  of  Eq. (2.4) 
is  satisfied  also.   In  order  to  select  the  physically 
realizable  solution  among  the  weak  solutions,  some 
additional  condition  has  to  be  satisfied.   The  criterion 
used  is  that  the  'correct'  solution  should  be  the  one 
obtained  in  the  limit  as  e    -*    0  in  the  viscous  equation 


(2.5)         ut  +  f(u)x  =   eu 


XX' 
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Lax  [28]  has  defined  an  entropy  inequality,  with  the  aid  of 
an  entropy  function  U(u)  for  Eg. (2.1),  which  is  reguired  to 
satisfy  the  following  two  conditions: 
(i)  U  satisfies 

(2.6)  Uu  fu  =  Fu. 

where  F  is  some  function  called  entropy  flux, 
(ii)  U  is  a  convex  function  of  u.  i.e., 

(2.7)  Uuu  >  0. 

Multiplying  Eg. (2.5)  by  Uu,  using  the  conditions  of 

Eg.  (2.6)  and  Eg.  (2.7),  also  taking  the  limit  as  e    ■+    0  shows 

that  the  correct  solution  is  the  one  satisfying 

(2.8)  U(u)t  +  F(u)x  <  0. 

If  u  is  a  weak  solution  with  discontinuity  travelling  with 
speed  s,  then 

(2.9a)        s« [U]  <  [F] . 

or 

(2.9b)       s. (UR  -  UL)  <  (FR  -  FL). 
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When  applied  to  the  gas  dynamics,  Eq.  (2.8)  or  (2.9)  is 
equivalent  to  the  fact  that  entropy  increases  across  the 
shock  waves.   For  this  reason  we  shall  call  Eq.(2.8)  or 
Eq. (2.9)  the  entropy  condition  or  entropy  inequality. 
In  the  scalar  case,  Oleinik  [29]  gave  the  entropy  condition 
across  the  shock  as 


f(u)  -  f(uL)          f(u)  -  f(uR) 
(2.10)         <  s  <  


u  -  uL  u  -  uR 

for  all  u  between  ul  and  ur. 

Notice  that  the  role  of  the  entropy  condition  is  to 
attempt  to  distinguish  those  discontinuous  solutions  which 
are  physically  realizable  from  those  which  are  not. 

Lax  [28]  then  defined  the  kth  field  to  be  genuinely 
nonlinear  if 

(2.11)  (grad  /cjj.ri  *  0. 

where  ri  is  the  kth  right  eigenvector  of  A(u) .   He  also 
defined  the  field  to  be  linearly  degenerate  if 

(2.12)  (grad  «i).ri  =  0. 

For  a  genuinely  nonlinear  field,  we  call  the  discontinuity 
a  shock  or  shock  wave  if  for  some  index  i 
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(2.13)  «i(uR)  <  s  <  /cj.(uL)  . 

For  linearly  degenerate  fields,  we  call  the  discontinuity  a 
contact  discontinuity  if 

(2.14)  *i(uR)  =  s  =  /ci(uL). 

In  one  spatial  dimension,  the  inviscid  equations  of 
gas  dynamics  can  be  written  in  the  conservative  form  as 


dq               3f(q) 
+  =  0, 

at     ax 
where 

q  =  (p,  m,  e  )T. 


is  the  vector  of  conservative  variables, 

f  =  (m,  (m2/P)+p,    (e+pjm/p  )T 

is  the  flux  vector,  and  m  =  pu.   In  the  above  equation  p  is 
the  pressure,  P    is  the  density,  u  is  the  velocity,  and  the 
total  energy  per  unit  volume  is  denoted  by  e  =  p/(7~l)  + 

m 

The  eigenvalues  of  the  Jacobian  matrix  A(q)  =  af/3q 
are 
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«i  =  u  -  c,    «2  =  U,     K3  =  u  +  c 

where  c  is  the  local  speed  of  sound  given  by  c  =  (7P/p)^. 
The  corresponding  eigenvectors  are 

ri  =  (1,  u-c,  H-uc)T 

r2  =  (1,  u,  hn2)T 

rj   =  (1,  u+c,  H+uc)T 

where  H  =  (e+p  )/P   =   c2/(7-l)+  hu2   is  the  enthalpy.   The 

wave  families  associated  with  k±   and  k3  can  be  rarefaction 

waves  or  shock  waves.  The  wave  associated  with  k2   is  a 

contact  discontinuity.   The  contact  discontinuity  is  due  to 

the  original  discontinuity  in  the  data;  it  is  presented  in 

the  linear  approximations  to  the  equations.   Therefore,  the 

corresponding  characteristic  field  is  called  linearly 

degenerate.   The  shock  waves  are  the  discontinuities  due  to 

the  nonlinearities  in  the  equations;  the  corresponding 

characteristic  field  is  called  genuinely  nonlinear. 

We  now  consider  an  explicit  (2k+l) -point  finite 

difference  scheme  in  conservation  form  approximating 

Eq. (2.1) 

n+1    n 
(2.15)      ui   =  ui  -  A[  hi+1/2  -  hi_1/2  ] 
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where 

A  =  At/Ax 

hi+l/2  =  h(ui-k+l/  ui-k+2/---#  ui+k) 
hi-l/2  =  h(ui-k/  ui_k+1,...,  Ui+k-x) 

and  hi+i/2,  the  numerical  flux,  must  be  consistent  with  the 
physical  flux  in  the  sense  that  h(u,u,...,u)  =  f (u) . 

In  1960,  Lax  and  Wendroff  [1]  showed  that  if  v(x,t)  is 
a  solution  of  such  a  finite  difference  scheme  in 
conservation  form,  and  if  v(x,t)  converges  almost 
everywhere  to  some  function  u(x,t)  as  ax  and  At  tend  to 
zero,  then  u(x,t)  is  a  weak  solution  of  Eq. (2.1). 

Conservation  forms  maintained  by  the  finite  difference 
schemes  are  very  important  since  this  ensures  that  the 
solutions  generated  by  such  schemes  have  the  correct  shock 
speed. 

Another  important  theoretical  contribution  was  made  by 
Harten  et  al.,  [6],  who  defined  a  class  of  monotone  finite 
difference  schemes.   A  scheme 


n+1 
(2.16)       ui    =  H(ui_k,  Ui_k+1,...,  Ui+k) 

n 
=  Ui  -A  [h(ui_k+1,. ..,  Ui+k) 

-  h(ui_k,  .  .  .,  ui+j^!)] 
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is  said  to  be  monotone  if  H  is  a  monotonically  increasing 
function  of  each  of  its  arguments.   All  convergent 
solutions  of  such  finite-difference  monotone  schemes  will 
satisfy  entropy  conditions  and  hence  are  physically 
realizable  solutions.   Harten  et  al.,  [6]  also  showed  that 
all  monotone  schemes  in  conservation  form  can  only  be 
first-order  accurate. 

2.2  Background  of  TVD  Schemes 

Consider  the  scalar  hyperbolic  conservation  law 

(2.17)  ut  +  [f(u)]x  =  ut  +  a(u)ux  =  0. 

A  weak  solution  to  this  problem  has  the  following 
monotonicity  properties: 

(1)  No  new  local  extrema  in  u(x)  may  be  created. 

(2)  The  value  of  a  local  minimum  is  nondecreasing  and 
the  value  of  a  local  maximum  is  nonincreasing. 

It  follows  from  these  monotonicity  properties  that  total 
variation  (TV)  in  x  of  u(x,t)  does  not  increase  in  t;  i.e., 

(2.18)  TV(u(t2)  )  <  TV  (u(t2)  ),    for  all  t2  >  tx. 

n 

Let  ui  be  the  numerical  approximation  of  the  solution  of 

Eg. (2.17)  at  x  =iAx  and  t  =  nAt,  with  Ax  the  spatial  mesh 
size  and  At  the  time  step.  Consider  a  one-parameter  family 
of  five-point  difference  schemes  in  conservation  form 
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n+1       n+1      n+1 

(2.19)  ui   +  \6(hi+1/2   -   hi_1/2) 

n  n        n 

=  ui  -  A(l-0) (hi+i/2  -  hi_1/2) 

n     n 
where  e    is  a  parameter,  A  =  At/Ax,  hi+1/2  =  h(ui_i  n\, 
n     n  ' 

ui+l/  ui+2  ) /  and  function  hi+i/2  is  called  the  numerical 
flux  function.  Let 

n         n+1 
Ri+l/2  ■  (l-0)hi+1/2  +  *hi+l/2 

be  another  numerical  flux  function.  Then  Eq.(2.19)  can  be 
rewritten  as 

n+1    n 

(2.20)  Ui   =  ui  -  A(fii+1/2  -  fii_1/2) 

This  new  numerical  flux  now  is  a  function  of  eight 
variables, 

B     n   n     n     n+1   n+1   n+1   n+1 
Ri+l/2  =  h(ui_!,  ui,  ui+1,  ui+2,  ui.i,  ui   ,  ui+1,  ui+2), 

and  is  consistent  with  the  physical  flux  f  in  the  following 
sense 

R(u,  u,  u,  u,  u,  u,  u,  u  )  =  f(u) 
To  simplify  the  notation,  Eq.(2.19)  can  be  rewritten  as 
(2.21a)      L.un+1  =  R.un 
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or 

(2.21b)      un+1  =  (L_1.R).un  =  S.un 

where  L  and  R  are  the  finite  difference  operators: 

(2.22a)      (L.u)i  =  ui  +  A0(hi+1/2  -  hi_1/2) 

(2.22b)      (R-u)i  =  ui  -  A(l-0) (hi+1/2  -hi_1/2) 

and 

S  =  L_1.R 

The  total  variation  of  a  mesh  function  un  is  defined  to  be 

(2.23)  TV(Un)  =  I       | ui+1  -  Ui|  =   |   |Ai+1/2un| 

i=— oo  i=— oo 

Here  and  throughout  this  studies,  the  general  convention 

Ai+l/2  (*)  =  (*)i+i  ~  (*)i 

is  used.  The  numerical  scheme  of  Eq. (2.19)  for  an  initial- 
value  problem  Eq. (2.17)  is  said  to  be  TVD  if 

(2.24)  TV(un+1)  =  TV(S«un)  <  TV(un) 

In  addition,  a  scheme  is  called  monotonicity 
preserving  if  the  finite  difference  operator  S  is 
monotonicity  preserving;  that  is,  if  un  is  a  monotone  mesh 
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function,  so  is  the  S»un.   Following  Harten  [11],  the  TVD 
schemes  are  monotonicity  preserving.   This  is  the  reason 
why  a  TVD  scheme  will  not  produce  spurious  oscillations. 

The  following  sufficient  conditions  for  Eq. (2.19)  to  be 
a  TVD  scheme  are  due  to  Harten  [11], 

(2.25a)      TV(R.un)  <  TV(un) 


and 


(2.25b)      TV(L.un+1)  >  TV(un+1) 

Assume  that  the  numerical  flux  h  in  Eq. (2.19)  is  Lipschitz 
continuous  and  that  Eq. (2.19)  can  be  written  as 

n+l    _-  + 

(2.26)    Ui   -A0(Ci+1/2Ai+1/2u  -  Ci-1/2Ai-1/2U)n+1 

n         .-  _  + 

-  Ui  +  A (1-0) (Ci+1/2Ai+1/2U  -  Ci_1/2Ai_1/2U)n 

where  Ci+1/2  =  c  (u^,  Ui,  ui+1,  ui+2)  and  Ci_1/2  =  C+( 
ui-2/  ui-if  ui,  Ui+1)  are  some  bounded  functions.  Then 
Harten  [11]  further  showed  that  sufficient  conditions  for 
Eq. (2.24)  are 
(a)  if  for  all  i 

±  ± 

(2.27a)      Ci+1/2  =  A(l-0)Ci+1/2  >  0 
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+  .+ 

(2.27b)   Ci+1/2   +  Ci+1/2  =  A (1-0) (Ci+1/2  +  Ci+1/2  )  <  1 

and 

(b)  if  for  all  i 

_  + 
(2.28)       -oo  <  c  <  -  X9    Ci+1/2  <  0 

for  some  finite  C.  This  set  of  sufficient  conditions  are 
very  useful  in  checking  or  constructing  the  TVD  schemes 
which  do  not  exhibit  the  spurious  oscillations  associated 
with  the  conventional  second-order  shock-capturing  schemes. 


CHAPTER  III 
UPWIND  AND  SYMMETRIC  TVD  SCHEMES 

In  this  chapter,  the  description  of  nonlinear,  upwind  and 
symmetric,  explicit  and  implicit,  second-order  accurate, 
entropy  satisfying  schemes  for  hyperbolic  conservation  laws 
will  be  presented.   The  notion  of  TVD  schemes  was  first 
introduced  by  Harten  [5].  He  developed  a  class  of  first-order 
TVD  schemes  that  prevent  the  growth  of  the  total  variation  of 
the  solution,  and  therefore  the  possibility  of  spurious 
oscillations  can  be  eliminated.   He  also  devised  second-order 
accurate  TVD  schemes  by  adding  limited  anti-diffusive  flux  to 
a  first-order  accurate  entropy  satisfying  upwind  scheme.   His 
method  is  sometimes  referred  to  as  modified  flux  approach.   A 
more  complete  description  of  the  TVD  scheme  based  on  Harten "s 
modified  flux  approach  will  be  presented  in  the  following 
sections. 

3.1  First-Order  Accurate  TVD  Schemp 

Consider  the  one-dimensional  scalar  hyperbolic 
conservation  law 
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au     af(u) 
(3.1)   +  
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au 

+   a(u) 

du 

•  —   0 

at      ax      at        ax 
and  with  the  initial  data 

(3.2)  u(x,0)  =  u0(x) 

where  a(u)  =  df/du  is  the  characteristic  speed. 

As  described  in  Chapter  II,  the  discrete  approximation  of 
Eq. (3.1)  by  a  three-point  explicit  conservative  finite 
difference  scheme  can  be  written  as 

n+1     n      n        n 
<3-3>      ui    =  »1  "  Mhi+i/2  "  hi_1/2) 

where  A  =  At/Ax. 

The  numerical  flux  in  Eq(3.3)  can  be  defined  as 

(3.4)     hi+1/2  =  ^[fi  +  fi+1  -  ^(ai+1/2)Ai+1/2u] 

where  fL  =   f(Ui)  and  Ai+1/2u  =  ui+1  -  Ui/  and 


(3.5)     ai+1/2  =  i 


(fi+l_fi)/(Ai+i/2u)     if  Ai+1/2u  *  0 


df/du  i  if  Ai+1/2u  =  o, 

U  =  Ui 


In  equation  (3.4),  rj,   can  be  viewed  as  the  coefficient  of 
numerical  viscosity,  and  it  is  a  function  of  A  and  ai+1/2. 
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Since  the  solution  of  the  difference  scheme  of  Eq.(3.3) 
does  not  necessarily  possess  the  TVD  property,  spurious 
oscillations  (wiggles)  may  be  produced,  especially  in  the 
neighborhood  of  discontinuities  (shock  waves,  contact 
discontinuities) .   In  order  to  eliminate  such  spurious 
oscillations  of  the  numerical  solution,  Harten  [5]  has  shown 
that  if  a  finite  difference  scheme  can  be  written  in  the  form 
of 

n+l    n     + 
(3.6a)     Ui    =  Ui  -  ACi_1/2Ai_1/2U  +  ACi+1/2Ai+1/2U 

where 

+ 
(3.6b)      ACi_1/2  =  A[  ai_1/2  +  tf(ai_1/2)]/2 

(3.6c)     ACi+1/2  =  a[  -ai+1/2  +  iKai+1/2)]/2 

then  the  sufficient  conditions  for  the  numerical  scheme  to  be 
TVD  are 


(3.7) 

ACi+l/2 

+ 

> 

0 

(3.8) 

ACi+l/2 

> 

0 

+ 

(3'9>  A(Ci+l/2  ♦  Ci+i/2)  *  l- 

Note  that  Eq. (3.9)  imposes  a  Courant-Friedrichs-Lewy  (CFL)- 
like  condition  on  the  difference  scheme  Eg. (3.3). 

Unlike  the  monotone  schemes,  TVD  schemes  do  not 
automatically  satisfy  the  entropy  condition,  and  the  schemes 
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may  converge  to  a  non-physical  (entropy  condition  violated) 
solution.   Consequently,  some  mechanism  may  have  to  be  added 
to  the  TVD  schemes  in  order  to  enforce  the  selection  of  the 
physical  solution.   As  shown  in  [11],  a  slight  modification  of 
the  coefficient  of  the  numerical  viscosity  term 


r  m 


(3.10) 


*(z)   =  i 


1*1  > 


(Z2  +  e2  )/2 


Z    <   £ 


can  avoid  violating  the  entropy  condition.   In  Eq. (3.10)  r/>(z) 
is  an  entropy  correction  to  |z|,  and  i  is  a  small  positive 
parameter  (see  reference  [11]  for  detailed  formulation  of  e). 

It  is  well  known  that  an  explicit  scheme  subject  to  a 
CFL-condition  is  limited  to  use  of  a  small  time  step. 
Therefore,  in  order  to  obtain  steady-state  solutions,  it 
usually  requires  long  computing  time.   This  penalty  is 
particulary  severe  for  high  Reynolds  number  Navier-Stokes 
calculations  where  the  fine  meshes  needed  to  resolve  small 
length  scales  presented  in  the  flow  field  dictates  the  use  of 
very  small  time  steps.   The  implicit  schemes,  on  the  other 
hand,  can  usually  adopt  much  larger  time  steps.   In  [18], 
Harten  devised  a  one-parameter  family  of  conservative  schemes 
of  the  form 


(3.11) 


n+1      n+1      n+l 
Ui  +  A0(hi+1/2  -  hi_1/2)  = 
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n         n        n 
ui  -A (1-0) (hi+1/2  -  hi_1/2) 


When  0=0,  Eq. (3.11)  is  reduced  to  Eq. (3.3),  an  explicit 
method.   When  6*0,    Eq. (3.11)  is  an  implicit  scheme.   For 
example,  if  6    =   1/2,  the  time  differencinq  is  the  trapezoidal 
formula,  and  if  9    =1,  the  time  differencinq  is  the  backward 
Euler  method.   In  qeneral  Eq.(3.11)  is  first-order  accurate  in 
space.   When  6    =  1/2,  the  scheme  is  second-order  accurate  in 
time.   Equation  (3.11)  can  be  written  as 

(3.12)  L«un+1  =  R»un 

where  L  and  R  are  the  finite-difference  operators,  defined  as 

(3.13a)         (L.u)i  =  ui  +  A0(hi+1/2  -  hi_1/2) 

(3.13b)         (R-u)i  =  ui  -  A(l-0)(hi+1/2  -  hi_1/2). 

A  sufficient  condition  for  Eq.(3.11)  to  be  a  TVD  scheme  is 
that 

(3.14a)  TV(R.un)  <  TV(un) 

(3.14b)  TV(L«un)  >  TV(un+1). 
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A  sufficient  condition  for  Eq. (3.14)    is  the  CFL-like  condition 


1 
(3.15)  Uai+l/2l    *   AV>(ai+1/2)    < 


i  -  e 

It  is  noted  that,  according  to  Eq.(3.15),  the  backward  Euler 
implicit  scheme  (0  =  1  in  Eq. (3.11))  is  unconditionally  TVD, 
while  the  trapezoidal  formula  (8   =  1/2)  is  TVD  under  the  CFL- 
like  restriction  of  2.   The  forward  Euler  explicit  scheme  {6    - 
0)  is  TVD  under  the  CFL  restriction  of  1. 

3.2  Second-Order  Accurate  TVD  Scheme 

In  order  to  achieve  a  second-order  accurate  scheme  which 
satisfies  the  TVD  conditions,  Harten  modified  the  numerical 
flux  fi  of  a  first-order  scheme  by  adding  a  limited  anti- 
diffusive  flux  gi  to  the  numerical  flux 

(3.16)  t±   -  t i .+  gi 

where  g^  approximates 


au 

1/2  Axy>(a)  

ax 


and  thus  cancels  the  first-order  error.   Therefore,  Harten 's 
method  is  sometimes  referred  to  as  the  modified-flux  approach. 
The  modified  flux  is  chosen  so  that  the  scheme  is  second-order 
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in  smooth  regions  and  first-order  at  points  of  extrema. 

After  introducing  the  modified  numerical  flux,  the  five- 
point,  second-order  accurate  backward  Euler  implicit  TVD 
scheme  (i.e.,  f  =  1)  of  Eq. (3.11)  can  be  written  as 

n+1    _n+l     _n+l       n 
(3.17)         ui  +  A(hi+1/2  -  hi_1/2)  m   ui 

and  the  modified  numerical  flux  hi+1/2  can  be  defined  as 
(3.18a)   hi+1/2  =  [fi  +  fi+1  -  \Hai+1/2  +  7i+l/2) M+l/2u]/2 


(gi+l-9i)/(Ai+i/2U)     if  Ai+1/2u  *  0 


(3.18b)  7i+l/2  ■ 


if  a 


1+1/2"  =  0 


where  £i  =  fi  +  gi 

It  is  necessary  to  limit  gi  to  prevent  the  possibility  of 
Ti+l/2  becoming  unbounded.   A  simple  choice  of  g^  is 

gi  =  S.max[0,  min(ai+1/2 |Ai+1/2u| ,  S.ai_1/2Ai_1/2u) ] 
where 

S  =  sign(Ai+1/2u) 
and 
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a  (Z)  =  **V>(Z)  >  0 

for  steady-state  calculation. 
One  can  rewrite  g^  as 

gi  =  minmod(ai+1/2Ai+1/2u,  <7i-i/2M-l/2u) 
where  the  minmod  function  is  defined  as  follows: 

minmod (x,y)  =  sgn(x) •max{0,  min[|x|,  y«sgn(x)]} 
or 


minmod (x,y)  =  • 


0  if  x  and  y  have  opposite  signs 

x  if  |x|  <  |y|  and  x  and  y  have  the  same  sign 

y  if  |y|  <  |x|  and  x  and  y  have  the  same  sign 


The  minmod  function  is  used  here  mainly  to  avoid  overshoots 
and  undershoots  in  the  vicinity  of  discontinuities.   Therefore 
gi  is  the  so-called  'flux  limiter'  function  for  the  control  of 
unwanted  oscillations  in  the  numerical  scheme. 

After  careful  examination  of  Harten's  modified-flux 
approach,  Yee  found  that  the  numerical  flux  hi+i/2  used  in  Eq. 
(3.18a)  is  quite  diffusive  [12].   A  slightly  modified  less 
diffusive  numerical  flux  can  be  defined  as  follows: 

(3.19)  hi+1/2  =   [fi  +  fi+1  +  *i+i/2] 
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where 

^i+1/2  -   <T(ai+l/2)  (9i+gi+l)-  ^(ai+i/2  +  7i+i/2)Ai+l/2u 
and 

a(z)  =  1/2* (z) 
gi  =  minmod(Ai+1/2u,  Ai_!/2u) 


7i+i/2  ■  CT(ai+l/2) 


(gi+l-gi)/(Ai+i/2U)  if  Ai+1/2  *  0 


if  Ai+l/2  =  0 


This  modified  form  is  that  the  ai+i/2  has  been  taken  out  from 
the  definition  of  g^. 

3.3  Extension  of  Scalar  TVD  Schemes  to  a  System  of 
Conservation  Laws 

There  are  two  difficulties  in  extending  the  scalar 
conservation  law  to  either  the  systems  of  eguations  or 
eguations  in  more  than  one  space  dimension.   First,  the  total 
variation  of  the  solution  to  the  system  of  hyperbolic 
equations  is  not  necessarily  a  monotonic  decreasing  function 
of  time.   The  total  variation  may  actually  increase  whenever 
the  waves  interact.   Second,  It  has  been  shown  by  Goodman  and 
Leveque  [30]  that  a  TVD  scheme  in  the  two  dimensional  case  is 
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at  most  first-order  accurate.   Neverthless,  It  is  possible  to 
extend  the  scalar  scheme  to  system  cases  such  that  the 
resulting  scheme  is  TVD  for  the  "  locally  frozen"  constant 
coefficient  system.   The  idea  behind  this  extension  relies 
heavily  on  the  Roe's  approximate  Riemann  solver  [30]. 
Consider  a  one  dimensional  hyperbolic  system 


au         3E(U) 
(3.20)         +   =  0 


at        ax 


or 


au         au 

(3.21)  +  A(U)  =  0 

at         ax 


Here  U  and  E(U)  are  column  vectors  of  m  components. 
Roe  first  linearized  Eq.(3.21)  to 


au    .   au 
+  a  =  o 


at       ax 

where  A(UR,UL)  is  a  constant  matrix  for  any  given  Ul,  Ur  which 
represent  the  local  conditions.   The  initial  data  for  this 
Riemann  problem  is 


U(x,0)  = 


'    UL    X  <  0 
.  UR    X  >  0 


The  matrix  A(UL,UR)  satisfies  the  following  properties  which 
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are  the  so-called  Property  U  defined  by  Roe  [15]. 

(1)  A(UL,UR)  is  a  linear  mapping  from  vector  space  U  to  the 
vector  space  E. 

(2)  As  UL  -  UR  +   U,  A(UL,UR)  *  A(U)  ,  where  A(U)  =  3E/3U. 

(3)  For  any  UL,  UR,  A(UL,UR) (UL-UR)  =  EL  -ER. 

(4)  The  eigenvectors  of  A  are  linearly  independent. 

The  eigenvalues  of  A  may  be  viewed  as  the  wave  speeds  of 
the  Riemann  problem  and  the  projection  of  (UL-UR)  onto  the 
eigenvectors  may  be  regarded  as  the  jump  between  the 
intermediate  states.   Conditions  (3)  ensure  that  the  exact 
solution  is  obtained  if  (UL,UR)  satisfy  the  jump  condition 

S(UL  -UR)  =  (EL  -ER) 

where  S  is  the  shock  speed.  Condition  (3)  shows  that  S  is  an 
eigenvalue  of  A.   Since  A  is  the  mean  value  Jacobian  ,  in  one 
dimensional  Euler  eguations  of  gas  dynamics,  such  a  matrix  is 
constructed  by  evaluating  its  Jacobian  matrix  with  specially 
averaged  cell  interface  values  (denoted  by  subscript  i+1/2)  of 
density,  velocity  and  total  enthalpy: 

Pi+1/2  =  (Pi)h(Pi+i)is 

ui+l^i+l)^  +  ui(pi)h 


ui+l/2  = 


(Pi+l)h   +  (Pi)* 
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Hi+l/2  = 


(Pi+l)h    +  (Pi)'5 


from  which  the  speed  of  sound  can  be  calculated  as 

2 

Ci+l/2  =  <(7-l)[  Hi+1/2  -  ^(ui+1/2)2]} 

Roe's  approximate  Riemann  solver  is  further  generalized  by  Yee 
[25].   Her  approach  is  sometimes  referred  to  as  the  "  local 
characteristic  approach  " .   The  procedure  is  to  define  at  each 
point  a  local  system  of  characteristic  variables  W  and  to 
obtain  a  system  of  uncoupled  scalar  equations.   After  that, 
one  can  apply  TVD  algorithm  of  a  scalar  hyperbolic 
conservation  law  to  each  of  the  locally  defined  (frozen 
coefficients)  characteristic  variables  of  Eq. (3.20)  as 
follows: 


(3.22)         +  A  =  0, 


aw 

+ 

A 

aw 

at 

ax 

w  = 

R' 

-1 

U 

A  =  R_1UR  =  diag(a1) 

where  diag(a1)  denotes  a  diagonal  matrix  with  diagonal 
elements  a1  which  are  the  eigenvalues  of  the  matrix  A(U) ,  R 
denotes  the  matrix  whose  columns  are  eigenvectors  of  A,  and 
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R-1  is  the  inverse  of  R. 

A  two-dimensional  hyperbolic  system  of  conservation  laws 
is  considered  here  to  demonstrate  the  construction  of  the 
scheme , 

(3.23)       Ut  +  E(U)X  +  F(U)y  =  0 

Here  U,  E(U)  and  F(U)  are  column  vectors  of  m  components.  Let 
A  =  aE/aU  and  B  =  aF/au  be  the  Jacobian  matrices.   Let 

the  eigenvalues  of  A  be  (ax,  ax, ,ax)  and  the 

eigenvalues  of  B  be  (ay,  ay, ,ay) .  Denote  Rx  and  Ry  as  the 

matrices  whose  columns  are  the  eigenvectors  of  A  and  B,  and 

-1     -1 
denote  Rx  and  Ry  as  the  inverses  of  Rx  and  Ry.   Let  the  grid 

spacing  be  denoted  as  Ax  and  Ay  such  that  xi  =  iAx  and  yj  = 

I)Ay.  Let  ai+1/2,  Ri+i/2>  Ri+l/2  denote  guantities  ax,  Rx,  Rx 

evaluted  at  Ui+1/2,j  respectively.   Similarly,  let  aj+1/2/ 

—1  p  _i 

Rj+l/2»  Rj+l/2  denote  the  guantities  ay,  Ry,  Ry  evaluated 

at  Ui^j+i/2  respectively. 

Define 


(3.24a)       ai+l/2  =  Ri+l/2  (ui+l,j  "  ui,j) 


as  the  difference  of  the  characteristic  variables  in  the  local 
x-direction,  and  define 


(3.24b)        "j+l/2  =  Rj+l/2  (ui,j+l  "  Ui,j) 
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as  the  difference  of  the  characteristic  variables  in  the  local 
y-direction. 

With  the  above  notation,  a  one-parameter  family  of  TVD 
schemes (3. 12)  in  two  dimensions  can  be  written  as 


n+1    x  _n+l       _n+l 
(3.25a)     U^j  +  A  0(Ei+1/2/j  -  Ei_1/2,j) 

y  _n+l       _n+l 
+  A  *(Fifj+i/2  "  Fi^j.i/2) 

n      x       _n         .n 
=  "i,j  -  A  (1  -0(Ei+i/2/j  -  Ei_1/2,j) 


y       .n         _n 
-  A  (1-  *)(Fifj+i/2  -  Fi^-i/2) 


with  Ax  =  At/Ax,  Ay  =  At/Ay.  Then  the  numerical  flux 
function  ^x+1/2,-}    for  hoth  the  upwind  and  symmetric  TVD 
schemes  can  be  expressed  as 

(3.25b)   Ei+1/2/j  =  %£«ifj  +  Ei+lfj  +  Ri+i/2*i+l/2]- 

Similarly,  one  can  define  the  numerical  flux  F^  j+i/2- 

3.4  Upwind  TVD  Schemes 

The  elements  of  $i+i/2/  needed  in  Eq.(3.19)  for  a  second- 
order  upwind  TVD  scheme,  originally  developed  by  Harten  and 

later  modified  and  generalized  by  Yee  [12],  are  denoted  by 
(^i+1/2)  1  £    "  1 , -  -  - ,m  and  defined  as  follows: 


(3.26a)  (    *i+i/2)U   = 
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I  I 

(ai+l/2)  (9i  +  <?i+l   ) 


i  i  i 

-  ^(ai+1/2   +  7i+i/2)Qi+l/2 


with 


(3.26b) 


a(z)    =   **(z) 


(4i+i  -  h 


2  a 

i)/ai+l/2     ai+l/2  *   ° 


(3.26c)    7i+l/2    =   CT(ai+l/2)    < 


ai+l/2 


=    0 


Examples  of  the  • limiter1  function  gi  can  be  expressed  as 


(3.27a) 


H  Z  i 

gi  =  minmod(ai+1/2/    «i-i/2) 


(3.27b)    gi  = 


it  s.  £ 

(ai+l/2ai-l/2    +    lai+l/2   ai-l/2l) 


ai+l/2   +  "i-1/2 


a  an 

(3.27c)  gi  =  S.max[    0,    min(2 |ai+1/2 | ,    S»ai_1/2) , 

J  i 

min( |ai+1/2| ,    2S«ai_1/2)] 


S  =  sgn    (ai+1/2) 
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Equation (3. 27b)  and  Eq. (3.27c)  are  suqgested  by  Roe[31]  and 
van  Leer [32]  respectively. 

3.5  Symmetric  TVD  Schemes 

In  1984,  Davis  [23]  showed  that  Roe-Sweby's  second-order 
TVD  scheme  can  be  interpreted  as  a  Lax-Wendroff  scheme  plus  an 
upwind  weighted  numerical  dissipation  term.  A  specific  flux 
limiter  then  can  be  constructed  by  satisfying  the  sufficient 
conditions  of  TVD  properties  mentioned  previously.   Also  the 
requirement  of  upwind  weighting  can  be  removed.   Thus  TVD 
schemes  are  not  necessarily  upwind  any  more,   these  non-upwind 
TVD  schemes  are  referred  to  as  symmetric  TVD  schemes.   Shortly 
after  that,  Roe  [24]  reformulated  Davis's  approach  in  such  a 
way  that  the  formulation  is  easier  to  analyze.  Flexible 
choices  of  flux  limiters  are  also  provided  to  enhance  the 
shock  resolution.   Subsequently  Yee  [25]  further  generalized 
the  works  of  both  Davis  and  Roe  to  a  class  of  symmetric  TVD 
schemes.   The  formulation  of  Roe-Davis  can  be  viewed  as  one  of 
the  variants  of  the  symmetric  TVD  schemes  proposed  by  Yee. 
The  advantages  of  Yee's  approach  are  that  (1)  stiff  problems 
can  be  treated  properly  by  using  implicit  schemes,  (2) 
numerical  boundary  conditions  are  not  so  sensitive  in 
comparison  with  upwind  TVD  schemes,  (3)  computational 
complexity  can  be  further  reduced. 
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I 

S 


The  elements  of  the  ($1+1/2)  denoted  by  (^i+1/2)   for  a 
general  second-order  symmetric  TVD  scheme  [25]  are 


I  i        i         J! 

(3.28a)       (^i+1/2)   ■  -V>(ai+l/2)  ["i+1/2  "   Qi+1/2  I« 


Examples  of  the  'limiter1  function  Qi+1/2  can  be  expressed  as 


i  1  J  I  J 

(3.28b)    Qi+1/2   =  minmod(ai_i/2,ai+i/2)    +  minmod(ai+i/2,ai+3/2) 

I 

"  ai+l/2 

I  2  I  I 

(3.28c)   Qi+1/2  =  minmod(ai_i/2,  ai+l/2f  ai+3/2) 

1  a  a  a 

(3.28d)      Qi+1/2   =  minmod    [    2ai_i/2,    2ai+1/2,    2ai+3/2, 

I  s. 

l/2(a-i-i/2    +   ai+3/2)] 

i  i 

where  ai+1/2  and  ^(ai+i/2)  are  defined  as  before.   In  general, 

the  limiter  of  Eq. (3.28d)  is  more  compressive  than  the  other 
two. 

Since  Eq(3.25a)  is  a  highly  nonlinear  implicit  scheme,  in 
order  to  gain  computational  efficiency,  local  linearizations 
are  applied  to  the  nonlinear  terms  (denoted  by  n+1  time  level) 
of  Eq. (3.25a).   There  are  two  linearized  forms  of  Eq. (3.25a) 
proposed  by  Yee.   The  first  method  will  maintain  the 


42 
conservative  property,  but  the  resulting  scheme  may  no  longer 
be  unconditionally  TVD.   This  method  is  referred  to  as  the 
linearized  conservative  implicit  (LCI)  form.   The  second 
method  will  destroy  the  conservative  property  but  preserve  the 
unconditionally  TVD  property.   The  latter  method  is  referred 
to  as  the  linearized  non-conservative  implicit  (LNI)  form. 
However,  according  to  Yee  and  Harten  [12],  It  appears  that  the 
LCI  form  had  better  convergence  rate  than  the  LNI  form.   Thus 
the  only  LCI  form  is  adopted  here. 

3.6  Linearized  Conservative  Implicit  (LCI)  Form 

The  nonlinear  terms  of  Eg. (3.25a)  are  linearized  in  time 
about  Un  by  a  Taylor  series  such  that 

En+1  =  En  +  An(un+1  _  ^  +  0(At2). 

Next  we  can  locally  linearize  the  coefficient  of  Ai+1/2  U  by 
dropping  the  time  level  from  (n+1)  to  n,  then  obtain  the  LCI 
form  in  delta  formulation  as 


XX  XX 

(3.29a)  [I   +   A    *Hi+1/2,j    -   A    *Hi_1/2,j 


y  y 

+    Ay^Hi/j  +  1/2    -    Ay<?Hifj_1/2](Un+l    -   U«) 


_n  _n 

=    -    Ax[Ei+1/2,j    -   Ei_1/2,j     ] 
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_n         _n 

^y[Fi,j+l/2  "  ^1^-1/2  ] 


where 


x  x 

(3.29b)      Hi+1/2/j  =  ^[Ai+1/j  +  Oi+l/2,j  ]n 


Y  Y 

(3.29C)      Hi/j+1/2  =  ^tBi,j+i  +  ni,j+l/2  ]n 


Here  A  and  B  are  the  Jacobians  of  the  fluxes  E  and  F,  and 


(3.30a)   Ai+i/2,j  =  (  Rx  diag[  p*-   iKa*  +  y*) ]Rx*) 1+1/2*1+1/2 

y  I    i   l   -1 

(3.30b)   ni^+i/2  =  (  Ry  diag[  0   -  *(a  +  7  )  ]Ry  )j+i/2Aj+i/2 


where 


(3.30c)       4+1/2  =  (si  +  gi+l)/«  1+1/2  » 

i        i,         2  a 

£j+l/2  =  (9j  +  gj+l)/"j+l/2 

For  steady-state  application,  one  way  to  simplify 
Eq.(3.30)  is  to  use  a  spatially  first-order  implicit  operator; 
e.g.,  by  redefining  Eq(3.30a)  and  Eq. (3.30b)  as 
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(3.31a)        «i+i/2,j   =    (Rx  diag[-  Vfa1)]    Rx1) i+i/2Ai+i/2 
(3.31b)        n?,j+i/2  =    (Ry  diag[-  ^(a1)]   Ry1) j+i/2Aj+i/2 


x 


The  computation  can  be  reduced  even  more   if  fii+i/2,j    and 

v  ... 

fli+j+1/2  are  simplified  to  diagonal  matrices.   For  example, 

one  can  redefine  Eg (3. 31)  as 


(3.32a)    0i+i/2,j  =  (  diag[  -  max  ^(a1)  ] ) i+i/2Ai+i/2 


y  . 

(3.32b)    "i,j+i/2  =  (  diag[  -  max  ^(a1)  ] ) 1+1/2* j+1/2* 

I 


Eguation(3.29)  together  with  Eg (3. 32)  is  referred  to  as  the 
linearized  conservative  diagonal  form. 

3.7  Alternating  Direction  Implicit  (API)  Form 

The  integration  of  the  two-dimensional  implicit 
operators,  such  as  Eg. (3.29a),  is  still  very  expensive.   One 
can  simplify  the  solution  process  by  introducing  an 
approximate  factorization  of  the  two-dimensional  operator  into 
two  one-dimensional  operators.   This  ADI  form  of  Eg. (3.29a) 
reads 
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X  X 

(3.33a)      [  I  +Ax0Hi+1/2,j  -  Ax0Hi_1/2,j  ]  D* 

=  -  Ax[Ei+1/2,j  -  Ei-i/2,j] 
"  AYtFi,j+i/2  "  Fi,j-l/2] 


y         y 

(3.33b)     [  I  +  Ay^Hi^+i/2  -  ^0^1,^-1/21    D  =  D* 


(3.33c)     Un+1  =  Un  +  D 

The  solution  algorithm  now  consists  of  two  one- 
dimensional  sweeps,  one  in  x-  and  one  in  y-direction.   Each 
sweep  requires  the  solution  of  a  linear  system  involving  a 
block  tridiagonal  matrix  which  can  be  solved  very  efficiently 
by  block  LUD  (lower-upper  decomposition)  solver. 


CHAPTER  IV 
FLUX  LIMITERS 

As  already  mentioned,  first-order  accurate  schemes,  in 
general,  suffer  from  numerical  diffusion,  white  the  classical 
second-order  accurate  schemes  such  as  the  Lax-Wendroff  scheme 
[1],  exhibit  the  spurious  oscillations  near  discontinuities  of 
the  solution.   In  order  to  overcome  this  problem,  One  can 
derive  the  flux  limiters  where  the  anti-diffusive  fluxes  are 
limited  so  as  to  prevent  these  oscillations  from  occurring. 
In  this  chapter,  the  mechanisms  of  a  class  of  flux  limiters 
are  analyzed.   Some  specific  flux  limiters  which  have  been 
proposed  by  various  authors  are  also  examined. 

Consider  the  initial  value  problem  for  a  scalar  hyperbolic 
conservation  law 


3u       3f(u) 
(4.1)       +  =  0, 

at      ax 


and  initial  value 

u(x,0)  =  u0(x) . 


As  before,  the  general  three-point  explicit  finite  difference 
in  conservation  form  which  approximates  Eq.(4.1)  can  be 
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written  as 


n+1     n      n        n 
(4.2)      ui    =  ui  -  A(hi+1/2  -  hi_1/2) 


where 


satisfying 


hi+l/2  =  h(ui+1,ui) , 


(4.3)  h(u,u)  =  f(u), 

that  is,  the  numerical  flux  h(u,u)  should  consist  with  flux 

function  f (u) . 

If  the  scheme  of  Eq. (4.2)  can  be  written  in  the  form 

n+1    n  n  n 

(4.4)  Ui    =  Ui  -  Ci_1/2Ai_1/2U   +  Di+1/2Ai+1/2U  , 

where  Ci_i/2  and  Di+i/2  are  functions  of  un,  then  it  can  be 
shown  [5]  that  sufficient  conditions  for  the  scheme  to  be  TVD 
are  the  satisfaction  of  following  inequalities: 

(4.5)  Ci+1/2  >  0, 

(4.6)  Di+1/2  >  0, 

(4.7)  Ci+1/2  +  Di+1/2  <  1. 

For  easy  presentation,  the  following  scalar  linear 
equation  is  considered 
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3u          an 
(4.8)      +  a  =o       a  =  const.  >  0. 

at       ax 

In  [34],  Sweby  has  shown  how  the  second-order  Lax-Wendroff 
scheme  could  be  constructed  by  adding  an  anti-diffusive  flux 
to  the  first-order  upwind  scheme,  that  is 

n+1     n         n  n 

(4.9  )   ui   =  Ui  -  ^Ai_1/2u  -  a_[  H   v(l-v)Ai+1/2u   ], 

where  CFL  number  w   -  Aa,  and  a_  is  the  backward  difference 
operator,  i.e.,   a_u  =  ui  -  *£«.].. 

It  is  well  known  that  the  Lax-Wendroff  scheme  produces 
spurious  oscillations  near  the  discontinuities  of  the 
solution,  while  the  first-order  upwind  scheme  does  not.   By 
limiting  the  amount  of  the  anti-diffusive  flux  in  Eq. (4.9),  it 
is  possible  to  obtain  the  second-order  accurate  TVD  scheme. 
The  concept  of  this  anti-diffusive  flux  is  much  the  same  as 
the  flux  corrected  transport  (FCT)  approach  proposed  by  Boris 
and  Book  [9],  and  Zalesak  [10].   The  main  differences  between 
the  present  flux  limiters  and  FCT  are  that  the  former  are 
single  step  schemes,  the  limiter  being  functions  of  data  un, 
whereas  FCT  is  a  two-step  procedure. 

The  anti-diffusive  flux  (the  second  term  on  the  right  hand 
side  of  Eg. (4.9))  is  modified  by  introducing  a  limiting 
function  <j>^, 
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(4.10)      -A_[  h    <t>i   v  (l-v)Ai+1/2U  ] 

The  ^i  is  chosen  to  be  a  function  of  the  ratio  of  consecutive 
cell  gradients  of 


Ai-1/2U 
(4.11)     a  =   tf(ri)   where  ri  = 


Ai+1/2U 

This  parameter  r^  plays  a  central  role  in  most  TVD  schemes, 
The  resulting  scheme  of  Eq.(4.9)  can  be  rewritten  as 


n+1    n  n 

(4.12)  Ui  =  Ui  -  u{l   +  %(l-v)[*(ri)/ri  -  4>{ri-i)-\)Li-1/2u   • 


This  is  a  scheme  in  the  form  of  Eg. (4.4)  with 

(4.13)   Ci_1/2  =  *<l  +  h(Z-»){4(v±y/ri  -   ^(ri_i)  ]  }, 
Di+l/2  =  0. 

Sufficient  conditions  for  the  scheme  of  Eg. (4.12)  to  be  TVD 
are  that  v   <  1  and 

(4.15)       U(ri)/ri  -  #(rj_i)|  <  2. 

Moreover,  Sweby  specified  that  <f>  (r)  >  0  for  r  >  0  and  that 
(f>  (r)  =0  for  r  <  0.   Under  these  additional  restrictions,  the 
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bound  of  Eq.    (4.15)    becomes 

(4.16)  0  <  <t>{r)/r  <   2;  0  <  <j>  (r)    <   2. 

Note  that  if  <f>  (r)  =  1,  Eq.(4.12)  reduces  to  the  central 
difference  Lax-Wendroff  scheme,  and  if  4>  (r)  =  r,  Eq.  (4.12) 
reduces  to  the  second-order  upwind  scheme  proposed  by  Warming 
and  Beam  [35] . 

The  region  defined  by  Eq. (4.16)  is  shown  in  Figure  4.1. 
Another  constraint  for  the  scheme  to  be  second-order  accurate 
is  imposed  by  requiring  that  the  amount  of  limited  flux  added 
to  the  first  upwind  scheme  be  a  weighted  averaged  of  the 
second-order  centered  Lax-Wendroff  scheme  and  the  second-order 
upwind  scheme  of  Warming  and  Beam.   The  corresponding  region 
with  this  extra  constraint  is  shown  in  Figure  4.2. 

In  [32],  van  Leer  averaged  two  second-order  schemes  to 
produce  a  second-order  TVD  scheme.   He  combined  a  monotonicity 
preserving  but  non-conservative  version  of  the  Lax-Wendroff 
scheme  and  the  second-order  upwing  scheme  of  Warming  and  Beam 
to  give  a  conservative  monotonicity  preserving  version  of 
Fromm's  scheme  [33].   The  parameter  he  used  as  a  "  smoothness 
monitor  "   is  defined  by 
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Warming-Beam 
r)=r 


Figure  4.1  TVD  region 
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Figure  4.2  Second  order  TVD  region 
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*i+l/2 


Ai+1/2U 


Ai_ 


i-1/2 


u 


which  is  the  reciprocal  of  the  ri-   After  reformulating  the 
van  Leer's  scheme,  the  equivalent  flux  limiter  is 


(4.17a) 


*(r)   = 


r  +  r 


|r|  +  1 


2r 


(4.17b) 


+(r}  = 


r  >  0 


1  +  r 


r  <   0 


which  is  a  smooth  curve  and  lies  in  the  second-order  TVD 
region  of  Figure  4.3. 

In  [17],  Roe  proposed  a  second-order  accurate  monotonicity 
preserving  scheme  using  an  incremental  (or  fluctuating) 
approach.   This  involved  the  allocation  of  cell  flux 
difference  to  update  the  values  of  u  at  the  ends  of  the  cells 
to  obtain  a  first-order  upwind  scheme  at  the  (n+l)At  time 
level.  Then  a  fraction  of  the  increment  is  transferred  across 
the  cell  against  the  direction  of  flow  to  achieve  second 
order. 
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0  i  2 

Figure  4.3  Van  Leer's  limiter  Eq. (4.17) 
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The  resulting  scheme  can  be  proved  to  be  monotonic  under 
certain  restrictions  on  <f>  (r)  .   The  most  commonly  used  limiter 
is  the  "  minmod"  limiter  which  corresponds  to  the  lower 
boundary  of  the  TVD  region  in  Figure  4.4,  and  which  may  be 
written  as 


(4.18a) 


4>{r)    =  max  [  0,  min  (l,r)  ], 


or 


(4.18b) 


*(r)  = 


0         r  <  0 
[     min(l,  r)   r  >  0 


Furthermore,  Roe  [31]  has  experimented  with  various 
transfer  functions  and  proposed  his  more  compressive 
"superbee"  transfer  function,  which  when  written  as  a  limiter 
is  expressed  as 


(4.19a) 
or 


<f>(r)    =  max  [0,  min(2r,l),  min(2,r)], 


(4.19b) 


*(r)  P"« 


min(2r,l) 
min(2,r) 


r  <  0 

0  <  r  <  1 

1  <  r 


which  corresponds  to  the  upper  boundary  of  TVD  region  (see 
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Figure  4.4  Minmod  liniter,  Eq. (4.18) 
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Figure  4.5) 

In  general,  Roe's  "superbee"  limiter  is  the  most 
compressive  among  the  above  three  <$>  (r)  functions.   Limiter 
(4.17b),  due  to  van  Leer,  is  found  to  produce  a  slightly 
better  shock  resolution  than  minmod  limiter  function 
Eg. (4.18a).   m  certain  cases  [31],  "superbee"  limiter  is  too 
compressive  and  might  produce  a  non-physical  solution. 
Therefore  using  different  limiters  for  different 
characteristic  fields  could  be  a  good  compromise.   For  the 
one-dimensional  Euler  eguation  of  gas  dynamics,  the 
characteristic  fields  consist  of  two  nonlinear  fields  u±c  and 
one  linear  field  u.   One  can  further  improve  the  sharpness  of 
the  contact  discontinuities  by  using  a  slight  diffusive 
limiter  such  as  Eq. (4.17)  for  nonlinear  fields  and  using  a 
more  compressive  limiter  such  as  Eq. (4.19)  for  linear  fields. 

The  above  idea  can  be  extended  to  a  nonlinear  scalar 
equation  as  well  as  to  a  multi-dimensional  system  of 
conservation  laws. 

In  1984,  Davis  [23]  interpreted  the  upwind  TVD  finite 
difference  schemes  analyzed  by  Sweby  [34]  as  Lax-Wendroff 
schemes  plus  an  upwind  artificial  dissipation  term.   Since  the 
determination  of  upwind  directions  for  hyperbolic  systems  is  a 
complicated  procedure,  he  modified  this  artificial  dissipation 
term  in  such  a  way  to  make  this  term  independent 
of  the  wave  directions.   The  result  is  a  dissipation  term 
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4>(r) 


3- 


2- 


Roe's  superbee  limiter 


0  1  2  3 

Figure  4.5  Roe's  superbee  limiter,  Eq(4.19) 


- 
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which  is  based  on  TVD  constraints  and  which  does  not  contain 
any  upwind  weighting.   Shortly  thereafter,  Davis •s  work  was 
reformulated  by  Roe  in  a  way  which  is  easier  to  analyze.   The 
analysis  exhibited  a  class  of  TVD  schemes  which  was  not 
observed  by  Davis.   Roe's  scheme  can  be  written  as 

n+1    n 
(4.20)  Ui    =  Ui  -  hv  (l+j/)Ai_1/2U  -hv(l-v)Ai+1/2u 

-  %|*j  (1-1* I)  (i-Qi-i/2)Ai-i/2u 

-  h\u\  (1-M)  (l-Qi+l/2)Ai+l/2U. 

Here  the  first  line  represents  the  original  Lax-Wendrof f 

scheme,  and  the  other  terms  represent  conservative  dissipation 

terms.   The  function  Qi+1/2  depended  on  three  consecutive 

gradients  Ai_;iy2u,  Ai+1/2u  and  Ai+3/2u  is  of  the  form 

+ 
(4.21a)    Qi+i/2  =  Q(  ri+1/2,  ri+1/2), 


where 


7       Ai-V2u     +        Ai+3/2u 
(4.21b)   ri+1/2  =  ;  ri+1/2  - 


Ai+1/2U  Ai+1/2U 

Roe  further  showed  that  Q  can  be  chosen  in  such  a  way  to 
ensure  that  Eg. (4.20)  is  TVD,  That  is  Eg. (4.20)  can  be  of  the 
form  of  Eg (4. 7)  with 
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(4.22)  Ci_1/2    =   i/[l-^(l-i/)Qi_i/2    +^(l-^)Qi+l/2/ri+l/2    1# 
Di+l/2    =    0. 

if  a  >  0,  and  the  case  for  a  <  0  follows  in  a  similar  way, 
i.e. , 

(4.23)  Ci_1/2  =  0, 

Di+l/2  -  1*1  [1  -  %.(1  -  IH)  Qi+1/2 

+  H(l  -    M)Qi-l/2/*i-l/2  ]• 

If  one  assumes  that  both  Q  and  Q/r  are  always  positive,  then  a 
set  of  sufficient  conditions  for  Eq.  (4.20)  to  be  TVD  are 


2 
(4.24a)       Qi+l/2  < 


1  -\u\ 

2 


(4.24b)      Qi+1/2/  ri+1/2  < 


M 


+  2 

(4.24C)       Qi+1/2/  ri+1/2  <  - 


M 


Based  on  these  bounds,  some  flux  limiters  Q  can  be  defined  as 


(4.25a)    Q(r  ,  r+)  =  minmod(l,  r")  +  minmod(l,r+)  -1, 
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(4.25b)     Q(r~,  r+)  =  minmod(l,  r~,  r+) , 

(4.25c)    Q(r_,  r+)  =  minmod  [2,  2r",  2r+,  h  (r~  +  r+) ] . 

The  graphical  representations  of  these  three  limiters  for 
symmetric  TVD  schemes  are  shown  in  Figure  4.6(a,b,c). 
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o              r"! 

1 

l 
r+-l 

r"+r+-l 

r+ 

-1 

1 
r"-l 

0 

(a)  Q(r~,r  )  =  minmod(l,r~)  +  minmod(l,r  )  -  1 


■i 


(b)  Q(r~,r+) 


rainmod(l,r~,r  ) 


0  0 

(c)  Q(r  ,r  )  =  minmod(2,2r~,2r+,l/2(r~+r+)) 

Figure  4.6  Graphical  representation  of.   three  limiters, 
Eq(4.6  a,b,c) 


CHAPTER  V 

ON  THREE-DIMENSIONAL  SYMMETRIC  AND  UPWIND 

TVD  SCHEMES  IN  GENERALIZED  COORDINATES 

An  original  derivation  of  the  implicit,  second-order 
accurate  symmetric  and  upwind  TVD  schemes  for  three- 
dimensional  hyperbolic  systems  of  conservation  laws  in 
generalized  curvilinear  coordinates  will  be  presented  in 
detail  in  this  chapter. 

5.1  The  3-D  Euler  Equation  in  Generalized  Coordinates 

The  inviscid  three-dimensional  conservation  equations  of 
mass,  momentum,  and  energy  can  be  represented  in  a  flux- 
vector  form  that  is  convenient  for  numerical  simulations  as 


(5.1a) 

aq 

-  + 

3E(q) 

-  + 

3F(q) 

ay 

+ 

3G(q) 

dz 

at 

ax 

-  o, 

where 
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(5.1b) 


q  = 


p 

pU. 
pW 
pW 

e 


E  = 


p\l 

pu2  +  p 

puv 

,      F  = 

pUW 

u(e+p) 

pv 

pVU 

pv2   +   p 

pVW 

v(e+p) 


G  = 


p\l 

pWU 

pWV 

p\t2  +   p 
w(e+p) 


The  primitive  variables  of  Eq. (5.1b)  are  the  density  p,    the 
velocity  components  u,  v  and  w,  and  the  pressure  p.   The  total 
energy  per  unit  volume  e  is  related  to  p  by  the  equation  of 
state  for  a  perfect  gas 


(5.1c) 


p  =  (7-1) [e  -  hp  (u2  +  v2  +  w2)] 


where  7  is  the  ratio  of  specific  heats. 

A  generalized  coordinate  transformation  of  the  form 

(5.2)    r  =  t, 

f  =  £ (t,  x,  y,  z)  is  the  logitudinal  coordinate 

v   =  ri(t,  x,  y,  z)  is  the  circumferential  coordinate 

f  ■'■»  t  (t,  x,  y,  z)  is  the  near  normal  coordinate 

which  maintains  the  strong  conservation  law  form  of  Eq. (5.1a) 
is  expressed  as 


(5.3) 


aq 


dr 


+ 


at(g) 

3£ 


*F(q) 


dr, 


+ 


aft(q) 


ar 


=  o, 


where 


A 

q 
I 

I 

g 
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q/J, 

(It  q  +  Cx  E  +  ^y  F  +  |z  G)/J, 

(»?t  q  +  »/x  E  +  »?y  F  +  f|  G)/J, 

(ft  q  +  fx  e  +  ry  f  +  fz  g)/j, 


and 


3(x,  y,  z) 

J  =  ,  the  Jacobian  of  transformation. 

dtt,    r#  D 


The  strong  conservative  form  of  the  transformed  equations  is 
maintained  in  order  to  admit  shock  capturing. 

The  contravariant  velocity  components  along  £-,  r/-,    and 
f -coordinates  are  defined  as  follows 


(5.4a) 


U=?t+^U+^yV+^W, 


V  =  r/t  +  »?X  u  +  »?V  V  +  ^Z  w/ 


w  =  ft  +  Tx  u  +  fv  v  +  fz  w- 


Here  E,  F  and  G  of  Eg. (5.3)  becomes 


(5.5) 


E  = 


pU 

pV 

pW 

1 

puU+   £xp 

1 

pUV+»?xp 

1 

puW+fxp 

J 

pVU+    £yP 

pWU+    £zp 

(e+p)U-|tp 

,    F   = 

J 

pW+r/yP 
pWV+r/zP 

(e+p)V-»7tp 

,G   =   

J 

pVW+fyP 

pwW+rzp 
(e+p)W-ftp 

66 


Let  A  =  3E/aq,  B  =  8F/dq   and  C  =  dG/dq;    then  the  Jacobians  A, 
B  and  C   of  ft,  £  and  6  can  be  written  as 

(5.6a)  A  =  £t  i  +  £x  a  +  £y  b  +  £z  c, 
(5.6b)  B  =  vt  I  +  ^x  A  +  r)Y  B  +  nz  C, 
(5.6c)       C   =   ft  I  +  rx  A  +  fy  B  +  fz  C. 

The  detailed  derivation  of  three-dimensional  compressible 
inviscid  gas  dynamic  equationss  in  generalized  curvilinear 
coordinates  can  be  found  in  Pulliam  and  Steger  [36]. 

Let  c  =  (tP/p)^  be  the  local  speed  of  sound,  The 
eigenvalues  of  A  are 

(5.7a)     (a^,  a^,  ac ,  a^ ,  a?  )  =  (u,  U,  U,  U+k^c,  U-k^c) 
where  k^  =  (£x2  +  £y2  +  ^2^.  The  eigenvalues  of  B  are 

(5.7b)     (ar,,    ar,*    a% ,    an'    a?)  =  (V,  V,  V,  V+k^c,  V-k^c) 
where  k„  =  (r,x2   +  „y2  +  r,z2)^.      The  eigenvalues  of  t   are 

(5.7c)     (ar,  ar,  ar  ,  a*,  af  )  =  (w,  W,  W,  W+kfc,  W-kfc) 
where  kr  =  (rj£2  +  fy2  +  tm*)h m 
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12    3    4    5 
Furthermore,  let  R^  =  (R^,  R^  ,  R^  ,  R^  ,  R^ ) , 

R„  =  (rJ,  R2,,  rJ,  rJ,  rJ),  and  Rf  =  (rJ,  R2,  r£,  R*,  R^) 

be  the  matrices  whose  columns  are  eigenvectors  of  k,    &  and  6. 

—1   —1      —1 
Let  R^  7  R^ ,  and  Rf   be  the  inverses  of  R^ ,  R^ ,  and  Rf .   The 

forms  of  the  matrices,  RK   with  k   =  £,  r; ,  or  f   and  their 

inverse  can  be  written  as 

(5.8) 


K   - 


*x 

ky 

kz 

1 

1 

kxu 

ky\l-kzp 

kzU+kyp 

u+kxc 

u-kxc 

k*v+kzp 

/CyV 

kzV-kxp 

V+kyC 

V-kyC 

kxw-kyp 

kyW+kxp 

iczw 

W+kzC 

w-kzc 

kxq2/2+ 

P  (kZV-icyW) 

kyq2/2+ 

p  (kxw-kzu) 

kzq2/2+ 

P  (icyu-kxv) 

H+C0 

E-c'e 

(5.9) 


-1 

R*  = 


Kx(1~lDl)~ 

«Xb2u 

Kz/p  + 

~Ky/p+ 

-«xb2 

(«ZV-KyW)/p 

"Xb2v 

«xb2w 

/Cy(l-bl)- 

~Kz/p  + 

Kyb2V 

*x/p  + 

-Kyb2 

(«;xW-/czU)/p 

/Cyb2U 

Kyb2W 

Kz(1_bl)" 

Ky/p  + 

-/cx/p  + 

/czb2w 

-«zb2 

(»CyU-/CXV)/p 

*Zb2u 

«zb2v 

(b1-s/c)/2 

(kx/C- 

(Ky/C- 

(/cz/C- 

b2/2 

b2u)/2 

b2v)/2 

b2w)/2 

(b1+fl/c)/2 

-(/CX/C+ 

-(/Cy/C+ 

-(/cz/c+ 

b2/2 

b2u)/2 

b2v)/2 

b2w)/2 
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where 

"x  =  Kx/(Kx2   +  «y2  +  «z2)^/ 

Ky  =  «y/(«x2  +  (Cy2  +  kz2)'2  , 

«z  =  «z/(/cx2  +  *y2  +  «z2)^/ 

5     =  KXU  +  /CyV  +  /CZW, 

q2  =  u2  +  v2  +  w2^ 

H   =  c2/(7-D  +  q2/2, 
bl  =  *>2  (q2/2), 
t>2  =  (7-l)/c, 
C  =  (7P/p)*, 

Let  the  grid  spacings  be  denoted  by  A£  ,  Arj    and  Af  such 
that  $   =  iA£ ,  r/  =  jA»?  and  f  =  kAf  .  Denote  qi+i/2,j,k  as  some 
symmetric  average  of  qi+i,-),*  and  qi,j,k  (for  example, 
Roe's  average).   Let  a{+1/2,    Ri+i/2/  Ri+l/2  denote  the 
quantities  of  a^  ,  R^  ,  R?   related  to  A  evaluated  at  qi+i/2,j,k 

I  -1 

Similarly,  let  aj+1/2/  Rj+l/2f  Rj+l/2  denote  the  quantities 

of  a^  ,  R,j ,  R^  related  to  I  evaluated  at  qi,j+i/2,k  and  ak+l/2' 

Rk+l/2/  Rk+l/2  denote  the  quantities  of  a^ ,  Rr ,  rJT1  related 
to  £   evaluated  at  qi,j,k+l/2-   0ne  can  define 
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(5.10a)   ai+i/2  = 


Ri+l/2  (<3i+l,j,k  -  gi,j,k) 


^(Ji+l,j,k  ♦  Ji,j,k) 


as  the  component  of  (qi+i,j,k  "  <3i,  j  ,k)  (omitting  the  j  and  k 

indices)  in  the  locally  i-th  characteristic  £ -direction. 

Denote 

-1 
1                     Rj+l/2  (qi,j+l,k  "  <3i,j,k) 
(5.10b)   aj+1/2  -  *-= • — , 

*(Ji,j+i,k  -  Ji,j,k) 

as  the  component  of  (qi,j+i,k  "  3i,j,k)  (omitting  the  i  and  k 
indices)  in  the  locally  i-th  characteristic  r? -direction. 
And  denote 

tn  «**,   J      Rk+i/2  (qi,j,k+i  -  qjyjfk) 

(5.10c)   «k+l/2  = • 

*(Ji,j,k+l  -  Ji,j,k) 

as  the  component  of  (qi,j,k+l  -  qi,j,k)  (omitting  the  i  and  j 
indices)  in  the  locally  i-th  characteristic  f -direction. 

I 

The  explicit  expression   for  vector  <*i+i/2   of  Eq.  (5.10a) 

is  given  by 


(5.11) 
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1 
"i+1/2 


<*i+l/2 


<*i+l/2 


4 
ai+l/2 


5 
ai+l/2 


-/CyAi+1/2W+kzAi+1/2V+ 

2  .  2 

[ci+l/2Ai+l/2^_Ai+l/2P]«x/ci+l/2 


*xAi+l/2w-*zAi+l/2u+ 

2  .  2 

[ci+l/2Ai+l/2P~Ai+l/2P]«y/ci+l/2 


-"xAi+l/2v+KyAi+l/2u+ 

2  .  2 

[ci+l/2Ai+l/2^-Ai+l/2P]«z/ci+l/2 


{Pi+l/2ci+l/2(KxAi+l/2u+KyAi+i/2v+KzAi+l/2w)+ 

2  2 

[ci+l/2Ai+l/2P-Ai+l/2P] >/ci+l/2 

(/'i+l/2ci+l/2(*xAi+l/2u+«yAi+l/2v+«zAi+l/2w)  + 

2  2 

[ci+l/2Ai+l/2^_Ai+l/2P] )/ci+l/2 


where 


(5.12a) 


Ai+1/2P   ■  Pi+l,j,k  "  Pi,j,k, 


(5.12b) 


Ai+1/2U  =  ui+1,j,k  -  ui^^k, 


(5.12c) 


>i+l/2v  =  Vi+1/j/k  -  vi(j/kl 
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(5.12d)         Ai+1/2W  =  Wi+lfj/k  -  Wi^j^k, 

(5.12e)       M+1/2P  =  Pi+i,j,k  "  Pi,j,k* 

and  Roe's  form  of  the  averaging  in  the  £ -direction  is: 
(5.13a)       Pi+1/2     =   (*>i+i,j,k)*  (Pi,j,-k)Hi 


is  i^i     „.             X*  Ui+1'3'k  +  ?*>*}** 
(5.13b)      Ui+1/2,jfk  "  


with 


(5.13c)      vi+1/2,j/k  = 


(5.13d)      Wi+1/2/j/k  = 


X*  +  1 


X*  vi+1/j,k  +  vi/j/k 


X*  +  1 


x*  wj+l,j,k  +  wifjyk 
X*  +  1 


(5.13e)      Hi+1/2/j/k  =  


X*    +    1 


(5.13f)  Ci+1/2/j,k   =    (7    -    1)  [    Hi+1/2fj/k    -    h> 


(Ui+l,j/k  +  v?+lfj/k  +  wf+1/j/k)] 


(5.13g)  X*    =       (pi+i,j/k  /Pi,!,*)*, 

(5.13h)  H  =  yp/p(y-l)    +  h(u2   +  v2   +  w2) 
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It  is  noted  that  here  H  represents  the  total  enthalpy  and 
should  not  be  confused  with  the  numerical  flux  function  H. 

Therefore,  to  use  Roe's  averaging  for  the  £- 
differencing,  all  one  has  to  do  is  to  compute  n±+i/2   j  k» 
vi+l/2,j,k/  wi+l/2,j,k  and  ci+l/2,j,k  in  Egs.(5.7),  (5.8), 
(5.9),  (5.10)  and  (5.11).  Similarly  Roe's  averaging  can  be 
applied  for  the  »?-  and  f -directions. 

5.2  Algorithm  of  3-D  TVD  schemes  in  generalized  Coordinates 
With  the  above  notation,  a  one-parameter  family  of  TVD 
schemes  can  be  written  as 

(5.14)      qi/j,k   +  *    * [Ei+i/2,j,k  ~  Ei_1/2jfk] 


+  A^[Fn^+1/2/k  -  Fn^.1/2,k] 


+  ^n~Gi+,],-k+l/2    -   Gn^,k-l/2] 


q?,j,k  +    (1   "   'H«!+l/2,j,lC  "  En-i/2/j/k3 


+    (1    -    0[Fn/j+V2/k   -    ~F^tj-1/2l-k] 


+    (1    -    9)[GitjfK+1/2    -    Gifj/k_1/2] 
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where  S  is  a  parameter,  A*  =  Ar/A^,  \n=  At/Atj  and  A^=  Ar/Af.  A 
particular  form  of  the  numerical  flux  function  Ei+1/2  j  k  f°r 
both  the  upwind  and  symmetric  TVD  schemes  can  be  expressed  as 

(5.15)  Ei+1/2,j,k  =  %(  &it1tlt  *   ^i+l,j,k  ♦  Ri+1/2  *i+l/2]- 

Similarly,  one  can  define  the  numerical  fluxes  F^  j+1/2  k  and 
Gi,j,k+l/2- 

Upwind  TVD  Schemes 

I 

The  elements  of  the  $1+1/2  denoted  by  (^i+i/2)U  for  a 

second-order  upwind  TVD  scheme,  originally  developed  by  Harten 
[5],  and  later  modified  and  generalized  by  Yee  [12],  are 


I  I       J      J 

(5.16  )  (^i+i/2)u  =  ^(ai+1/2) (gi+i  -  gi) 

a  i.  1 

-  V>(ai+i/2  +  7i+i/2)<*i+l/2- 

For  steady-state  calculation  the  function  a(z)  =  l/2V»(z),  and 
the  functions  ip   and  7  are  defined  as  before,  i.e., 


(5.17)      rf>(z)    =  - 
and 


1*1  |z|  > 

.  (Z2  +  £2  )/2e  |z|  < 


e 


(5.18)7i+l/2  =  *(ai+l/2)  < 
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1  a       a 

(9i+l  "  gi)/«i+l/2 


"i+1/2  *° 


ai+l/2  =0 


Examples  of  the  "limiter"  functions  g^  can  be  expressed  as 

£  in. 

(5.19a)       gi  =  minmod  (a±-1/2i    «i+l/2) / 

a  a  sl  i.  i 

(5.19b)  gi    =    (ai+1/2    «i-l/2    +     l<*i+l/2    »i-l/2l)/ 

i  i 

(«i+l/2    +   ai-l/2)  / 

(5.19c)  gi  =  S.  max[0,  min(2  l«i+1/2 I ,  S'al-i/2), 


with 


min( |ai+1/2| ,  2S»ai_1/2)], 


S  =  sign  (ai+1/2) 


Again,  the  minmod  function  of  Eq. (5.19a)  is  defined  as  follows 


minmod(x,y)  =  sgn(x)«max{  0,min[  |x|,  y«sgn(x)]} 


Symmetric  TVD  Schemes 


The  elements  of  $1+1/2,  denoted  by  (^i+i/2)S  for  a 
general  second-order  symmetric  TVD  scheme  [25],  are 
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(5.20)  (^i+i/2)S  "  "  ^(ai+1/2)     t   «i+l/2   "  Qi+1/2    ] 


A 


Examples  of  the  "limiter"  functions  Qi+1/2  can  be  expressed  as 


(5.21a)       Qi+1/2  =  minmod(ai_1/2f  ai+l/2) 


+  minmod(ai+1/2,  <*i+i/2)  -  "i+1/2  > 


Ai  a  a  sl 

(5.21b)      Qi+1/2  =  minmod(ai_1/2,  <*i+i/2/  "i+1/2)  > 

k  4  J 

(5.21c)  Qi+1/2  -  minmod[2ai_1/2,    2ai+1/2/ 


I  i  i 

2ai+l/2    /    ^(ai-l/2   +  "i+1/2)] 


5.3  A  Conservative  Linearized  API  Form  for  Steady-State 
Application 

A  conservative  linearized  ADI  form  of  Eq. (5.14)   for 
steady-state  application  can  be  written  in  the  generalized 
coordinates  as 


i     i  i     i 

(5.22a)    [I  +  A  0    Hi+1/2,j,k  -  A  9    Hi_1/2,j/k]  D** 


i    ~n  _n 

"  A  [Ei+1/2/j,k  -  Ei_1/2/j/k] 
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~  A  [Fi,j+l/2,k  -  Fif j-i/2/k] 

f  -n  _n 

~  A  [Gi,j,k+l/2  ~  Gi/j/k_1/2]/ 


(5.22b)  [I  +  A  |  Hif j+i/2,k  "  A I  IlJ, j-l/2/k]  D*  =  D**, 


(5.22C)  [I  +  A  9    Hif  j/k-fl/2  "  A  I  Hif  j,k-i/2]  D  =  D*, 


(5.22d) 


qn+1  =  -n  +  D/ 


where 


(5.22e) 


Hi+i/2,j,k  =  *f*l+i,j,k  +  ni+i/2,j,k]n, 


9 

(5.22f)      Hi/j+1/2/k  =  %t«J,j+l,Jc  +  "i,j+l/2,k]n, 

r  r 

(5.22g)  Hi/j/k+1/2    =   ^[6i>j/k+1   +  Oi,j,k+1/2]n, 

and 

€  i 

(5.22h)      ni+1/2/j/k  =  diag[-  max  #(   ai+1/2)]    Ai+1/2, 


1  ^ 

(5.221)      ni,j+i/2/k  =  diagt-  max  iHaj+1/2)]    Aj+1/2, 

X* 

(5.22J)      ni,j/k+1/2   =  diag[-  max  tf(ak+1/2)]    Ak+1/2, 
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Here  &i+1/j/k,  &i,j+i,k>  and  ^i,j,k+l  are  Eq.(5.6)  evaluated 
at  (i+l,j,k),  (i,j+l,k)  and  (i,j,k+l)  respectively. 

Observe  that  Eq.(5.22)  is  the  original  Beam  and  Warming 

I  i  r 

algorithm  [3]  if  0i+i/2,j,k  =  "i,j+l/2,k  =  "i,j,k+l/2  =  0  in 
Eq. (5.22h)-(5.22j)  and  $i+i/2Ri+l/2  i-n  Eq(5.15)  is  replaced  by 
the  conventional  fourt-order  dissipation  term.   Therefore,  the 
implemention  of  this  conservative  linearized  ADI  scheme  into 
an  existing  central  difference  code  (such  as  the  one  based  on 
the  Beam  and  Warming  algorithm)  is  relatively  simple. 

One  has  only  to  add  the  extra  matrices  fli+i/2,j,k' 

n  f 

ni,j+l/2,k»  anc*  °i,j,k+l/2   for  tne  implicit  operator  and  uses 

a  more  sophisticated  dissipation  term  $i+i/2Ri+l/2  for  tne 
explicit  operator. 


CHAPTER  VI 
THIN  LAYER  NAVIER-STOKES  SOLVERS 

The  basic  governing  equations  and  numerical  algorithm 
used  in  this  chapter  have  been  adopted  from  the  works  of 
Steger  [37],  and  Pulliam  and  Steger  [36],   The  three- 
dimensional  Reynolds-averaged  Navier-Stokes  equations  in 
generalized  curvlinear  coordinates  are  simplified  by  using 
the  thin-layer  assumption.   These  equations  can  be  further 
reduced  to  a  set  of  axisymmetric  equations  which  govern 
azimuthal  invariant  flows.   The  resulting  r; -invariant 
equations  are  solved  by  using  an  implicit,  non-iterative 
approximate  factorization  finite  difference  scheme 
developed  by  Beam  and  Warming  [3].   The  axisymmtric  thin- 
layer  Navier-Stokes  code,  which  was  originally  developed  by 
the  U.S.  Army  Ballistic  Research  Laboratory  [26],  [38],  is 
modified  to  incorporate  the  TVD  schemes  derived  in  Chatpter 
V  and  to  adopt  a  spatially  variable  time  step  procedure  to 
enhance  both  the  numerical  accuracy  and  computational 
efficiency.   In  the  following,  a  brief  presentation  is  made 
to  describe  the  assumptions  and  simplifications  of  the 
governing  equations,  the  turbulent  closure  model  adopted, 
and  the  implementation  of  the  TVD  schemes. 
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6.1  Thin  Layer  Approximation 

In  external  high  Reynolds  number  viscous  flows,  the 
effects  of  viscosity  are  usually  concentrated  near  the  no 
slip  boundaries.  Hence  the  flow  profile  near  the  solid 
surface  must  be  accurately  resolved,  especially  in  the 
direction  normal  to  the  solid  surface.   In  order  to  provide 
an  adequate  spatial  resolution,  grid  lines  near  the  solid 
boundary  are  exponentially  clustered  (so  that  at  least  one 
or  two  grid  points  will  be  within  the  viscous  sublayer) . 
Without  flow  separation,  it  is  usually  assumed  that  only 
the  viscous  terms  in  the  ($•-)  direction  normal  to  the  solid 
surface  are  needed  to  be  retained  while  the  viscous  terms 
in  both  the  streamwise  (£-)  and  spanwise  (r,-)    directions 
can  be  neglected. 

The  resulting  equations  with  this  thin  layer 
approximation  can  be  written  as  follows  : 


aq      at  dP  a6     1    a§ 

(6.1)   +   +  +  = 


dT  d£  dr)  3f       Re    9i 

where  the  viscous  terms  in  the  f -direction  have  been 
collected  in  the  vector  |,   Detailed  derivation  for  both 
the  three  dimensional  full  Navier-Stokes  equations  and  the 
thin  layer  approximation  in  the  generalized  curvilinear 
coordinates  are  given  in  Appendix  A. 
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6-2  Azimuthal  for  n-)    Invariant  Equations 

Azimuthal-invariant  equations  are  obtained  from 
Eq. (6.1)  by  makinq  use  of  followinq  two  restrictions:  (l) 
all  body  qeometries  are  axisymmetric  to  the  r? -direction  and 
(2)  the  state  variables  and  the  contravariant  velocities  do 
not  vary  in  the  t? -direction.   The  equations  resultinq  from 
simplifications  were  derived  by  Nietubicz  et  al  [26],  [38] 
as  shown  in  the  followinq 


aq 

(6.2a)                      + 

3E 

a6 

l        a§ 

dr 

3£ 

a  r            Re      d  r 

where 

(6.2b) 
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(6.2c) 


ft   =   J-l 


and 


° 
pV[R£(U-£t)+Rr(W-ft)] 
-pVR(V-i,t)-p/R 

0 


*i  =  ^(rx2  +  fz2). 

m2  =  (m/3)  (rxur  +  rzwf ), 

m3  =    (u2   +  v2  +  w2)j./2   +  MPr~1(7-l)"1( 
u  =  ft  *  ^xu  +  ^zw, 
V  =  ft  +  ^yV/ 

w  "  ft  ♦  fxu  +  rzw- 


c   K' 


6.3  Turbulence  Model 

The  algebraic  mixing  length  model  developed  by  Baldwin 
and  Lomax  [39]  is  incorporated  into  the  present  thin  layer 
formulation  to  model  the  effect  of  turbulence.   This  is  a 
two-layer  algebraic  model.   The  inner  layer  is  governed  by 
the  Prandtl  mixing  length  with  the  Van  Driest  damping  term; 
the  outer  layer  follows  Clauser's  approximation.   Computed 
vorticity  is  used  in  defining  the  reference  mixing  length 
required  for  the  outer  layer.   The  present  turbulence  model 
is  appropriate  for  attached  and  mildly  separated  flow 
problems. 
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6.4  Implicit  Central  Difference  Algorithm 

The  numerical  scheme  used  here  is  an  implicit,  non- 
iterative  approximate  factored,  finite  difference  algorithm 
written  in  delta  form,  proposed  by  Beam  and  Warming  [3]. 
The  resulting  finite  difference  equations  can  be  expressed 
as 

n  n     n         n   n 

(6.3a)   L^Lr[  Aq  ]  =  -At(S^t   +    S^   -  Re_1«r§  +  ft  )  -  De, 

with  the  linear  operators  L^  and  Lr  as 

(6.3b)       L£  -  [I  +  hS^n  -  Dil^], 

(6.3c)       Lf  =  [I  +  h5rCn  -  hRe-^J"1  ftn  J  -  Dilr]. 

Here  h  =  At,  the  S's  are  typically  three  point  second- 
order  accurate  central  difference  operators  and  A  and  v  are 
forward  and  backward-difference  operators  respectively.   J 
is  the  Jacobian  of  the  coordinate  transformation.   The 
Jacobian  matrices  A,  C  and  ft  result  from  the  local  time 
linearization  and  can  be  found  in  Appendix  A.   In  the 
current  study,  the  first-order  accuracy  in  time  and  second- 
order  central  difference  in  space  are  applied  throughout 
the  code 

Since  central  discretized  factorization  scheme  was 
used,  the  high  frequency  error  caused  by  the  nonlinear 
terms  may  not  be  completely  damped  out,  which  can  lead  to 
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numerical  solutions  with  spurious  oscillations.   Therefore 
numerical  dissipation  terms  denoted  as  Di  and  De  have  been 
inserted  into  Eq. (6.3).   The  explicit  fourth-order 
artificial  dissipation  term 

(6.4)     De  =  ee  At  J-i^A^V^)2  +  (A^Vj-)2]  J, 

is  added  to  the  right-hand  side  of  Eq. (6.3)  and  the 
implicit  second-order  smoothing  terms 

(6.5a)     Di  |?  =  ei  At  J-1  (AV)^  J> 

(6.5b)     Di  |r  =  €±    At  J"1  (AV)r  J, 

are  inserted  into  the  respective  implicit  block  operators. 
The  implicit  second-order  smoothing  terms  are  chosen  to 
keep  the  left  hand  side  factors  of  the  block  tridiagonal 
structure.   The  parameters  e±   and  ee  are  chosen  such  that 
*i    *  2ee. 

The  solution  algorithm  now  consists  of  two  one- 
dimensional  sweeps,  one  in  the  f-  and  the  other  in  the  $■- 
direction.   Each  step  requires  the  solution  of  a  linear 
system  involving  a  block  tridiagonal  matrix,  which  is 
solved  by  a  block  LUD  (lower-upper  decomposition)  solver. 
The  resulting  solution  process  is  much  more  economical  than 
the  unfactored  algorithm  in  computer  storage  and  CPU  time. 
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Also  steady  state  calculation  can  be  achieved  by 
integrating  the  governing  equations  from  an  arbitrary 
initial  condition  to  a  time  asymptotic  state. 

6.5  Dissipation  Model  of  TVD  Scheme 

As  mentioned  in  Chapter  V,  the  major  difference 
between  the  Beam-Warming  scheme  and  the  TVD  scheme  is  that 
the  TVD  scheme  requires  more  sophisticated  dissipation 
terms.  The  dissipation  terms  of  the  Beam-warming  scheme  in 
Eq.(6.4)  and  Eq. (6.5)  are  replaced  by  the  following 
expressions 

(6.6)     De  =  -H    At(Ri+1/2$i+1/2  -  Ri-i/2*i-i/2 

+  Rk+l/2*k+l/2  "  Rk-l/2*k-l/2)/ 
and 

(6.7a)    Di  \(   =  -h   At  (ni+1/2/k  -  fli_1/2/k), 
(6.7b)     Di  | r  =  -h   At  (ni/k+i/2  -  ni,k-i/2K 


where  o 1+1/2 fk  and  ni,k+l/2  are  defined  similarly  as  in 
Eq. (5.22h)  and  Eq. (5,22i)  except  that  the  j  index  is 
ommitted. 
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6.6  Grid  System 

A  schematic  illustration  of  a  Secant-Ogive-Cylinder- 
Boattail  (SOCBT)  projectile  shape  is  shown  in  Figure  6.1. 
Specifically,  the  projectile  under  study  is  3-caliber  for 
secant-ogive  part,  2 -caliber  for  the  cylinder  part  and  1- 
caliber  7  degree  for  boattail  part  which  is  further 
extended  for  another  1.77  caliber  to  meet  straight  sting. 

A  grid  generation  code  utilizing  a  hyperbolic 
differential  equation  [40]  is  employed  to  generate  a 
network  of  constant  £  and  f  lines  in  the  physical  x-z  plane 
as  illustrated  in  Figure  6.2.   Corresponding  uniform  grids 
of  £  and  f  in  the  computational  space  define  a  one  to  one 
mapping  between  points  j,k  in  the  physical  plane  to  points 
j,k  in  the  computational  plane  also  shown  in  Figure  6.2. 
The  grid  spacing  in  the  computational  domain  is  chosen  to 
be  unity,  that  is  A£  =  l  and  Af  =1. 

Once  the  projectile  geometry  is  determined,  the  grid 
points  along  the  body  surface  can  be  redistributed  by  using 
a  cubic  polynomial  as  the  clustering  function  to  make  much 
denser  grids  near  the  high  gradient  region. 
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Pysical  Domain 
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Figure  6.2  Mapping  from  physical  domain  to  computational 
domain. 
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For  the  projectile  problem,  since  the  outer  boundary 
is  unconstrained,  it  only  needs  to  be  sufficiently  far  away 
from  the  inner  boundary.   Such  a  grid  generation  problem  is 
frequently  encountered  in  external  aerodynamic  application. 
In  the  current  study,  the  outer  boundary  is  taken  as  4 
projectile  lengths  away  from  the  projectile  body.   The  "  C- 
type  "  grids  are  then  generated  by  the  hyperbolic 
formulation  [41] .   When  solving  the  viscous  flow  problem  it 
is  necessary  to  have  a  very  fine  mesh  normal  to  the 
projectile  body  near  the  wall,  so  that  the  velocity 
gradients  within  the  boundary  layer  can  be  adequately 
resolved.   The  grid  spacing  can  be  smoothly  controlled  by 
using  an  exponential  clustering  function  in  the  z- 
direction.   The  minimum  grid  spacing  used  near  the  wall  in 
the  code  is  2xl0~5  diameters  of  projectile. 

6.7  Numerical  Boundary  Conditions 

The  general  curvilinear  coordinate  transformations  are 
made  such  that  physical  boundaries  are  mapped  to  boundaries 
in  the  computational  domain.   This  makes  the  formulation 
and  implementation  of  boundary  conditions  easier.   For  the 
axisymmetric  projectile  flow  problem  with  an  attached 
sting,  unknown  values  of  q  on  the  boundaries  are  updated 
explicitly  after  each  integration  step.   The  numerical 
boundary  conditions  for  the  projectile  calculations  are  as 
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follows: 

(1)  Along  the  outer  boundary,  line  AB  in  Figure  6.2,  free 
stream  values  are  specified. 

(2)  Along  the  line  AD  in  Figure  6.2,  the  flow  variables  are 
obtained  by  a  second-order  extrapolation.   Numerically  this 
is  done  by  using  a  three-point  formula  to  approximate  the 
first  derivative  with  respect  to  £ . 

(3)  Along  the  outflow  boundary,  line  BC  in  Figure  6.2,  all 
the  flow  variables  are  obtained  by  a  two-point 
extrapolation.   However,  for  the  subsonic  cases  the 
pressure  is  specified  along  the  outflow  boundary- 

(4)  Along  the  inner  boundary,  line  CD  in  Figure  6.2,  the 
no-slip  condition  is  employed  for  three  velocity  components 
(i.e.,  u,v,  and  w  =  0) .   The  surface  pressure  is  obtained 
from  the  normal  momentun  equation  3P/dn  =  0,  where  n  is  the 
local  normal  to  the  solid  surface.  The  surface  density  is 
obtained  from  an  adiabatic  wall  assumption.   Finally,  total 
energy  is  computed  from  the  equation  of  state. 

6.8  Spatial-Varying  Time  Step 

For  steady  state  flow  problems,  it  is  desirable  to 
achieve  the  converged  solution  as  quickly  as  possibble. 
The  convergence  to  a  steady-state  solution  can  be 
accelerated  by  using  spatially  varying  time  steps  [4].   For 
a  fixed  time  step,  the  Courant  number  is  not  uniform  since 
the  grid  spacings  vary  from  very  fine  to  very  coarse  in  the 
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whole  flow  field.   The  use  of  a  spatial  varying  time  step 
can  thus  be  interpreted  as  an  attempt  to  use  a  more  uniform 
Courant  number  throughout  the  field.   The  choice  of  this 
spatial  varying  time  step  is  governed  by  the  formula 


At  ref 
(6.8)  At  =  1 


1  +  (J)^ 
where  the  J  is  the  Jacobian  of  transformation. 


CHAPTER  VII 
RESULTS  AND  DISCUSSION 

7.1  Numerical  Implementation 

One  objective  of  the  present  study  is  to  compare  the 
conventional  Beam-Warming  scheme  with  the  TVD  scheme  using 
different  flux  limiters  in  terms  of  accuracy,  stability  and 
efficiency.   In  this  chapter  the  turbulent  transonic  flow 
field  about  a  secant-ogive-cylinder-boattailed  (SOCBT) 
projectile  model  is  adopted  as  the  test  problem.   An 
axisymmetric  thin-layer  Navier-Stokes  code  developed  at  the 
U.S.  Army  Ballistic  Research  Laboratory  [26,38]  has  been 
modified  in  the  present  study.   One  major  difference  between 
the  original  and  modified  codes  is  the  treatment  of  the 
numerical  dissipation  terms.   This  approach  of  implementing 
the  TVD  scheme  is  suggested  by  Yee  [42],  who  interpreted  the 
TVD  algorithms  as  a  form  of  nonlinear  numerical  dissipation. 
This  viewpoint  can  greatly  reduce  the  work  of  incorporating 
the  TVD  schemes  into  the  existing  computer  code.   In  the 
original  version  of  the  code,  as  reported  in  [26],  the  Beam- 
Warming  scheme  was  adopted.   Based  on  this  framework,  the  TVD 
schemes  can  be  incorporated  into  the  original  computer  code  by 
adding  a  more  sophisticated  numerical  dissipation  model.   With 
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this  improved  numerical  dissipation  mechanism,  one  can  strive 
both  to  eliminate  the  non-physical  oscillations  commonly 
encountered  in  the  numerical  solution  yielded  by  the  Beam- 
Warming  scheme  and  to  satisfactorily  resolve  the  sharp 
gradients  present  in  the  flow  field.   In  the  present 
investigation,  three  variants  of  the  flux  limiters  have  been 
assessed  in  detail;  two  of  them  are  designed  for  the  upwind 
TVD  scheme,  and  the  third  one  is  for  the  symmetric  TVD  scheme. 
Specifically,  the  flux  limiter,  Eg(5.19b),  of  the  upwind  TVD 
scheme  suggested  by  van  Leer  is  called  'Limitl';  the  one  which 
combines  Roe's  superbee  limiter,  Eq. (5.19c),  (in  the  linear 
field)  and  van  Leer's  limiter  ,Eq. (5.19b),  (in  the  nonlinear 
field)  is  called  'Limit2',  The  flux  limiter  of  the  symmetric 
TVD  scheme,  Eq. (5.21c),  is  called  'Limit3'. 

Both  the  original  and  modified  thin-layer  Navier-Stokes 
codes  have  been  used  to  compute  transonic  turbulent  flow  over 
the  SOCBT  projectile  model  at  zero  angle  of  attack.   The  test 
cases  considered  here  are  the  SOCBT  projectile  with  (a)  Me  = 
0.96,  (b)  Mo,  =  0.91,  and  (c)  M*,  =  1.2.   Here  M*,  is  the  free 
stream  Mach  number.   Figure  7.1  shows  a  90x60  "  C-type"  grid 
for  a  flow  domain  of  4  projectile  lengths  with  special 
clustering  on  both  the  ogive-cylinder  and  cylinder-boattail 
junctures.   Moreover,  an  exponential  clustering  function  is 
employed  such  that  the  minimum  grid  spacing  chosen  above  the 
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Figure  7.1  The  90x60  "C-type"  grid  for  the  SOCBT  projectile 
with  sting. 
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Figure  7.2a  The  surface  pressure  coefficients  of 
Beam-Warming  scheme  for  the  SOCBT 
projectile  M„  =  0.96,  and  Re  =  7.7xl05. 
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projectile  surface  in  the  f -direction  is  2xl0~5  diameters  of 
the  projectile  in  order  to  adequately  resolve  the  viscous 
sublayer.   For  the  original  thin-layer  code,  the  parameters  of 
user-controlled  smoothing  terms  are  set  to    e^  =   2ce  =  4At. 
For  the  modified  code,  the  coefficient  of  e  in  Eq.(3.11)  is 
set  to  zero;  also  the  spatially  varying  time  step,  Eq. (6.7), 
can  be  used  to  further  enhance  the  convergence  rate.   No 
special  effort  was  devoted  to  vectorizing  the  program  for 
either  version  of  the  code.   Hence  the  CPU  time  presented  for 
the  following  cases  should  be  viewed  only  as  reference.   All 
the  numerical  computations  were  run  on  the  CRAY-XMP/48 
supercomputer  at  the  National  Center  for  Supercomputing 
Applications  (NCSA) ,  at  the  University  of  Illinois  at  Urbana- 
Champaign. 

7.2  Results  of  Flow  of  H  =  0.96 

Figure  7.2a  shows  the  converged  solution  for  the  surface 
pressure  coefficient  Cp  calculated  by  the  original  Beam- 
Warming  scheme.   The  symbols  "  o  "  on  the  plot  indicate  the 
experimental  data  obtained  by  Kayser  et  al  [27],  while  the 
solid  line  represents  the  numerical  solution.   As  one  can  see 
in  Figure  7.2a,  the  numerical  result  compares  favorably  with 
the  experimental  data  except  for  the  peak  point  located  in  the 
shock  region  above  the  middle  of  cylinder  segment. 

The  corresponding  comparision  for  the  results  obtained 
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with  the  TVD  schemes  are  presented  in  Figure  7.2(b,c,d). 
Overall,  the  surface  pressure  coefficient  yielded  by  both  the 
Beam-Warming  and  TVD  schemes  are  of  similar  quality  on  the 
present  90x60  grid  system.   A  comparison  of  the  pressure  and 
shear  drag  coefficients  predicated  by  the  four  schemes  is  made 
in  Table  7.1.   The  nondimensionalized  pressure  drag 
coefficients,  defined  as 

Dp  =  (pressure  force) /(pM  a*,2  L2) 

obtained  by  all  four  schemes  are  very  close,  with  less  than  1 
percent  difference  between  any  two  predictions.   The 
nondimensionalized  shear  drag  coefficient,  defined  as 

Df  =  (shear  force)/ (0.5  Pa>   V^2) 

shows  more  discrepancies  among  the  predictions.   The  Beam- 
Warming  scheme  predicts  the  smallest  value  of  Df,  while  the 
symmetric  TVD  scheme  predicts  the  largest  shear  drag  among  all 
the  schemes  under  study.   The  difference  between  the  largest 
and  the  smallest  predicted  values  of  Df  is  more  than  two 
percent.   Although  there  are  no  experimental  measurements 
available  for  a  definite  assessment  with  respect  to  the  drag 
coefficients,  the  results  shown  in  Table  7.1  broadly  confirm 
the  effects  of  the  TVD  schemes  on  the  numerical  solution. 
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Figure  7.2b 


The  surface  pressure  coefficients  of 
upwind  TVD  scheme  with  Limit 1  for  the 
projectile  with  M*,  =  0.96,  and 
Re  =  7.7xl05  using  a  90x60  C  grid. 
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Figure  7.2c 


The  surface  pressure  coefficients  of 
upwind  TVD  scheme  with  Limit2  for  the 
projectile  with  Ma,  =  0.96,  and 
Re  =  7.7xl05  using  a  90x60  C  grid. 
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Figure  7 . 2d  The  surface  pressure  coefficients  of 
symmetric  TVD  scheme  with  Limit3  for 
the  projectile  with  Me  =  0.96,  and 
Re  =  7.7x10s  using  a  90x60  C  grid. 


98 


In  terms  of  the  computational  efficiency,  Figure  7.3 
shows  the  convergence  history  of  the  four  schemes  under  study. 
The  symmetric  TVD  scheme  has  the  fastest  convergence  rate. 
Furthermore,  while  there  exist  differences  in  the  rate  of 
convergence  among  the  three  TVD  schemes,  Figure  7.3  clearly 
indicates  that  their  perforamces  are  all  much  superior  to  that 
of  the  Beam-Warming  scheme.   The  other  advantage  of  the  TVD 
schemes  is  the  fact  that  the  user  does  not  have  to  resort  to 
empirical  rules  to  choose  the  damping  parameters.   In  the 
Beam-Warming  scheme,  the  proper  choices  of  the  free  parameters 
such  as  q  and  ee  are  not  always  clear;  a  trial-and-error 
approach  is  usually  needed.   In  this  regard,  the  TVD  schemes 
can  be  viewed  as  a  more  rigorous  and  sophisticated  approach  in 
determing  the  desirable  numerical  dissipation. 

Table  7.2  presents  more  information  of  the  computational 
aspects  of  the  four  schemes.  Only  fixed  time  step  is  used  by 
the  Beam-Warming  scheme,  whereas  TVD  schemes  are  used  both 
with  and  without  the  provision  of  the  spatially  varying  time 
step.   On  the  basis  of  the  CPU  time  per  iteration,  the  Beam- 
Warming  scheme  is  the  most  efficient.   This  is  clearly 
expected  since  the  TVD  schemes  need  to  make  a  substantial 
amount  of  extra  computations  to  eliminate  the  spurious 
oscillations  and  to  produce  high  profile  resolutions.   Among 
the  TVD  schemes,  the  symmetric  TVD  scheme  has  the  best 
performance  in  terms  of  CPU  time  per  iteration.   Combined  with 
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Figure  7.3   Convergence  history  of  different  schemes 

(a)  Beam-Warming  scheme 

(b)  Upwind  TVD  scheme  with  Limit2 

(c)  Upwind  TVD  scheme  with  Limit 1 

(d)  Symmetrtic  TVD  scheme  with  Limit3 
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the  superior  convergence  rate,  the  symmetric  TVD  scheme  is  the 
most  efficient  scheme  overall.   But  nevertheless,  Table  7.2 
demonstrates  that  with  the  identical  grid  distribution,  all  of 
the  TVD  schemes  can  yield  the  steady  state  solution  with  much 
less  CPU  time  than  the  Beam-Warming  scheme.   With  the  aid  of 
the  spatially  variable  time  steps,  significant  improvements  in 
convergence  rate  have  been  obtained.   The  symmetric  TVD  scheme 
is  the  most  attractive  one  for  the  present  turbulent  transonic 
flow  calculations. 

Table  7.1   Pressure  and  Shear  Drags  for  Different  Schemes  of 
90x60  Grid  at  *!«,  =  0.96 


Scheme 

Limiter 

Pressure 
Drag  x  10~3 

Shear 
Drag  x 

10"3 

B-W 

- 

40.96 

27.88 

UPTVD 

Limit 1 

40.56 

28.22 

UPTVD 

Limit2 

40.67 

28.14 

SYMYVD 

Limit3 

40.72 

28.52 
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Table  7.2a  L2-Norm  Residual  and  CPU  Time  for  Different 
Schemes  with  90x60  Grid  at  Mc  =  0.96  Using 
Fixed  Time  Step 


Scheme 

Limiter 

L2-Norm 
Residual 

NO.  Of 

Iter. 

CPU  per 
Iter. 

Total  CPU 
(min) 

B-W1 

- 

9.6xl0~5 

6000 

0.73 

72.14 

UPTVD2 

Limit 1 

5.8X10"5 

1700 

0.92 

25.00 

UPTVD 

Limit2 

7.1X10"5 

1700 

0.96 

27.00 

SYMTVD3 

Limit3 

7.2X10"5 

1300 

0.85 

17.65 

1  :  Beam-Warming  Scheme 

2  :  Upwind  TVD  Scheme 

3  :  Symmetric  TVD  Scheme 

Table  7.2b  L2-Norm  Residual  and  CPU  Time  for  Different 
Schemes  with  90x60  Grid  at  M,  =  0.96  Using 
Spatially  Varying  Time  Step 


Scheme 

Limiter 

L2~Norm 
Residual 

No.  of 
Iter. 

CPU  per 
Iter. 

Total  CPU 
(min) 

UPTVD 

Limit 1 

6.3xl0~5 

1000 

0.90 

15.18 

UPTVD 

Limit2 

5.1X10"5 

1400 

0.95 

22.35 

SYMTVD 

Limit3 

7.2X10""5 

700 

0.81 

10.72 

7.3  Grid  Density  and  Numerical  Accuracy 


In  Figure  7.2,  satisfactory  and  comparable  numerical 
accuracies  among  the  four  schemes  were  obtained  on  a  90x60 
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grid  system.   In  order  to  assess  the  effectiveness  of  the  TVD 
schemes  and  the  relative  merits  between  them  and  Beam-Warming 
scheme,  a  series  of  computations  has  been  conducted  by  using  a 
progressively  coarser  grid  distribution.   First,  the  number  of 
grid  points  along  the  streamwise,  i.e.,  £-,  direction  is 
reduced  from  90  to  60.   That  is,  the  grid  system  now  employs 
60x60  nodes  instead  of  90x60  nodes.   Figure  7.4  presents  the 
comparisons  between  the  experimental  results  and  numerical 
ones  in  the  same  flow  condition  on  this  grid  system  among  the 
four  schemes.   Compared  to  Figure  7.2,  the  solutions  obtained 
by  the  TVD  schemes  on  both  the  90x60  and  60x60  grid  system  are 
essentially  undistinguishable  in  term  of  the  pressure 
coefficient.   The  Beam-Warming  scheme,  on  the  other  hand, 
shows  degradation  in  accuracy  on  the  coarser  grid  system, 
notably  in  the  shock  region. 

Next,  the  number  of  grid  points  along  the  transverse, 
i.e.,  f-,  direction  is  reduced  from  60  to  40.   Figure  7.5 
compares  the  prediction  and  measurements  of  the  pressure 
coefficient  on  a  90x40  grid  system  for  the  Beam-Warming  scheme 
and  symmetric  TVD  scheme.   Overall,  the  Beam-Warming  scheme 
fails  to  capture  the  shock  adequately  and  its  accuracy  is  much 
inferior  to  the  results  obtained  by  using  the  TVD  schemes. 
Again  all  three  variants  of  the  TVD  scheme  perform  comparably 
well,  although  only  one  of  them  is  presented  in  Figure  7.5. 
The  impact  of  the  reduction  of  grid 
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Figure  7.4a 


The  surface  pressure  coefficients  of 
Beam-Warming  scheme  for  projectile 
with  M<*>  -  0.96,  and  Re  =  7.7x10s 
using  a  60x60  C  grid.   The  L2~norm 
residual  of  5.16xl0~5  can  be  reached  in 
around  5000  time  steps. 
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Figure  7.4b 


The  surface  pressure  coefficients  of 
upwind  TVD  scheme  with  Limitl  for 
projectile  with  M»  =  0.96,  and 
Re  =  7.7xl05  using  a  60x60  C  grid. 
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The  L2-norm  residual  of  3.9xlO~D  can 
be  reached  in  around  1500  time  Steps. 
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Figure  7.4c 
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The  surface  pressure  coefficients  of 
upwind  TVD  scheme  with  Limit2  for 
projectile  with  M»  =  0.96,  and 
Re  =  7.7x10s  using  a  60x60  C  grid. 
The  Lo-norm  residual  of  8.5x10-°  can 
be  reached  in  around  1900  time  steps, 
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Figure  7.4d 


The  surface  pressure  coefficients  of 
symmetric  TVD  scheme  with  Limit3  for 
projectile  with  M»  =  0.96,  and 
Re  =  7.7xl05  using  a  60x60  C  grid. 
The  L2-norm  residual  of  3.4xl0-5  can 
be  reached  in  around  1100  time  steps. 
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Figure  7 . 5a 


The  surface  pressure  coefficients  of 
Beam-Warming  scheme  for  projectile 
with  M«  =  0.96,  and  Re  =  7.7xl05 
using  a  90x40  C  grid.   The  L2-norm 
residual  of  2.59xl0"6  can  be  reached 
in  around  3000  time  steps. 
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Figure  7 . 5b 
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The  surface  pressure  coefficients  of 

symmetric  TVD  scheme  with  Limit3  for 

projectile  with  M«  =  0.96.  and 

Re  =  7.7xl05  using  a  90x40  C  grid. 

The  L2-norm  residual  of  4.8xl0-r>  can  be 

reached  in  around  1100  time  steps. 
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size  on  the  numerical  accuracy  at  the  present  level  of  grid 
size  is  not  visible  for  the  TVD  schemes.   Further  reductions 
of  the  number  of  grid  points  have  also  been  performed. 
Figures  7 . 6  and  7 . 7  show  the  comparisons  of  numerical  results 
and  measured  data  of  the  three  TVD  schemes  on  the  80x40  and 
70x40  grid  system  respectively.   Again,  the  three  schemes  all 
appear  to  be  able  to  maintain  good  accuracy  in  both  cases.   In 
summary,  by  reducing  the  grid  size  from  90x60  to  70x40,  the 
TVD  schemes  are  shown  to  be  both  robust  and  accurate  enough  to 
yield  a  steady-state  solution. 

Table  7 . 3  presents  some  comparisons  of  the  performance 
among  the  TVD  schemes  on  the  70x40  grid  in  the  computational 
aspect.   Overall  the  symmetric  TVD  scheme  is  fairly  efficient 
in  term  of  iterations  and  CPU  time  per  iteration.   However  the 
convergence  rate  for  symmetric  TVD  scheme  has  slightly 
degraded  as  shown  in  Figure  7.8.   Table  7 . 4  compares  the 
predicted  pressure  and  shear  drags  on  the  70x40  grid  system. 
It  appears  that  results  yielded  by  upwind  TVD  with  Limitl  and 
the  symmetric  TVD  with  Limit3  closely  follow  each  other,  while 
that  by  the  upwind  TVD  scheme  with  Limit2  predicts  smaller 
drag  forces. 
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Figure  7 . 6a 


The  surface  pressure  coefficients  of 
upwind  TVD  scheme  with  Limitl  for 
projectile  with  M«  =  0.96,  and 
Re  =  7.7xl05  using  a  80x40  C  gri£. 
The  L2~norn  residual  of  7.68x10"°  can 
be  reached  in  around  1300  time  steps. 
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Figure  7 . 6b 


The  surface  pressure  coefficients  of 
upwind  TVD  scheme  with  Limit2  for 
projectile  with  M=°  =  0.96,  and 
Re  =  7.7x10s  using  a  80x40  C  grid. 
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The  L2~norm  residual  of  8.15x10"°  can 
be  reached  in  around  1500  time  steps. 
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Figure  7 . 6c 


The  surface  pressure  coefficients  of 
symmetric  TVD  scheme  with  Limit3  for 
projectile  with  M<*>  =  0.96,  and 
Re  =  7.7xl05  using  a  80x40  C  grid. 
The  L2~norm  residual  of  5.4xl0~5  can 
be  reached  in  around  1100  time  steps. 
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Figure  7.7a  The  surface  pressure  coefficient  of 
upwind  TVD  scheme  with  Limitl  for 
projectile  with  M^  ■  0.96,  and 
Re  =  7.7xl05  using  a  70x40  C  grid. 
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Figure  7.7b  The  surface  pressure  coefficient  of 
upwind  TVD  scheme  with  Limit2  for 
Projectile  with  Me  =  0.96,  and 
Re  =  7.7xl05  using  a  70x40  C  grid. 
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Figure  7.7c  The  surface  pressure  coefficient  of 
symmetric  TVD  scheme  with  Limit3  for 
projectile  with  Mo  —   0.96,  and 
Re  =  7.7xl05  using  a  70x40  C  grid. 
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Figure  7.8   The  convergence  history  of  different 
TVD  schemes  at  M*  =  0.96,  using  a 
70x40  grid. 

(a)  Upwind  TVD  scheme  with  Limit2 

(b)  Upwind  TVD  scheme  with  Limitl 

(c)  Symmetric  TVD  scheme  with  Limit 3 
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Table  7.3   L2-Norm  Residuals  and  CPU  Time  for  Different 
Schemes  of  70x40  Grid  at  M*,  =  0.96 


Scheme 

Limiter 

L2~Norm 
Residual 

No.  of 
Iter. 

CPU  per 
Iter. 

Total  CPU 
(min. ) 

UPTVFD 

Limitl 

4.2X10"6 

1500 

0.44 

10.90 

UPTVD 

Limit2 

1.4X10"5 

1900 

0.49 

17.27 

SYMYVD 

Limit3 

3.2X10"5 

1100 

0.42 

7.83       ! 

Table  7.4  Pressure  and  Shear  Drags  for  Different  Schemes 
for  70x40  Grid  at  M*,  =  0.96. 


Scheme 

Limiter 

Pressure 
Drag  xl0~3 

Shear 
Drag  xl0~3 

UPTVD 

Limitl 

43.42 

28.09 

UPTVD 

Limit2 

41.57 

27.68 

SYMTVD 

Limit3 

43.20 

28.05 

7.4  Results  of  Flow  of  M^  =  0.91 

In  this  case  the  transonic  speed  is  reduced  from  M^  = 
0.96  to  Me  =  0.91.   All  the  calculations  were  conducted  on  the 
70x40  grid  system.   The  computational  aspects  of  the  TVD 
schemes  are  described  in  Table  7.5.   The  computed  pressure  and 
shear  drags  are  shown  in  Table  7.6.   Figures  7.9(a,b,c)  show 
that  the  agreements  of  computed  surface  pressures  with 
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Figure  7 . 9a  The  surface  pressure  coefficients  of 
upwind  TVD  scheme  with  Limit 1  for 
projectile  with  M,,  =  0.91,  and 

Re  =  7.7xl05  using  a  70x40  C  grid. 
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Figure  7.9b  The  surface  pressure  coefficient  of 
upwind  TVD  scheme  with  Limit2  for 
projectile  with  M»  =  0.91,  and 
Re  =  7.7xl05  using  a  70x40  C  grid. 
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Figure  7.9c 


The  surface  pressure  coefficients  of 
symmetric  TVD  scheme  with  Limit3  for 
projectile  with  Me  =  0.91,  and 
Re  =  7.7xl05  using  a  70x40  C  grid. 
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experimental  data  are  fairly  good.   Figure  7.10  shows  the 
convergence  history  of  the  TVD  schemes  with  three  different 
flux  limiters.   The  convergence  rate  is  degraded  when  the 
upwind  TVD  scheme  with  Limit2  is  used.   Roe's  "  superbee  " 
limiter  used  in  the  linear  field  tends  to  produce  more 
compressive  profiles  than  other  limiters.   The  results 
obtained  here  suggested  that  the  convergence  rate  will  degrade 
when  too  compressive  limiter  is  used  for  the  flow  with  weak 
shocks. 


Table  7.5  I^-Norm  Residuals  and  CPU  Time  for  Different 
Schemes  of  70x40  Grid  at  M^,  =  0.91 


Schemes 

Limiter 

L2 -Norm 

No.  of 

CPU  per 

Total  CPU 

Residual 

Iter. 

Iter. 

(  Min.) 

UPTVD 

Limitl 

1.3xl0~5 

1900 

0.47 

14.75 

UPTVD 

Limit2 

3.4xl0~4 

1300 

0.49 

12.29 

SYMYVD 

Limit3 

5.1X10"6 

1500 

0.42 

10.61 

Table  7.6   Pressure  and  Shear  Drags  for  Different  Schemes 
of  70x40  Grid  at  M*,  =  0.91. 


Schemes 

Limiter 

Pressure 
Drag  x  10~3 

Shear 
Drag  x 

10"3 

UPTVD 

Limitl 

21.82 

25.25 

UPTVD 

Limit2 

18.80 

23.07 

SYMTVD 

Limit3 

20.97 

25.29 
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Figure  7.10   The  convergence  history  of  different 
TVD  schemes  at  M^,  ■  0.91  using  a 
70x4  0  grid. 

(a)  Upwind  TVD  scheme  with  Limit2 . 

(b)  Upwind  TVD  scheme  with  Limitl. 

(c)  Symmetric  TVD  scheme  with  Limit3 . 
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7.5  Result  of  Flow  of  1L  =  1.2 

In  this  case,  as  the  Mach  number  is  increased  from  M,*,  = 
0.96  to  Ma,  =  1-2,  computations  on  the  70x40  grid  system  have 
been  conducted  for  the  three  TVD  schemes.   Figures  7.1l(a,b,c) 
show  that  surface  pressure  coefficients  computed  by  using  the 
TVD  schemes  are  in  excellent  agreement  with  the  measured  data. 

All  three  TVD  schemes  appear  to  produce  almost  identical 
results.   Figure  7.12  shows  the  convergence  histories  of  the 
three  TVD  schemes.   The  convergence  rate  of  the  symmetric  TVD 
scheme  has  shown  difficulties  in  reducing  the  residual  to  a 
very  low  level  whereas  the  other  two  schemes  can  rapidly 
reduce  the  residuals  throughout  the  whole  calculation.   It  is 
noted  that  for  all  the  calculation  with  TVD  schemes,  the 
coefficient  of  entropy  correction  e    is  assigned  to  be  zero  in 
the  current  code.   Yee  has  pointed  out  that  the  smaller  the  € 
being  used,  the  slower  is  the  convergence  rate  of  the 
symmetric  TVD  scheme  [42].   Although  the  convergence  rate 
could  be  further  improved  if  a  positive  value  of  e  is  used,  it 
was  decided  not  to  pursue  this  in  the  present  work  since  one 
of  the  major  attractions  of  the  TVD  scheme  is  to  be  free  from 
arbitrarily  choosing  the  adjustable  parameter  values.   The 
pressure  drag  and  shear  drag  of  this  case  are  shown  in  Table 
7.7.   The  predictions  of  the  two  drag  coefficients  are  very 
close  among  the  three  schemes. 
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Figure  7 . 11a  The  surface  pressure  coefficients  of 
upwind  TVD  scheme  with  Limitl  for 
projectile  with  Ma  =  1.2,  and 
Re  ■  7.7xl05  using  a  70x40  C  grid. 
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Figure  7.11b  The  surface  pressure  coefficients  of 
upwind  TVD  scheme  with  Limit2  for 
projectile  with  !!«,  =  1.2,  and 
Re  =  7.7xl05  using  a  70x40  C  grid. 
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Figure  7.11c  The  surface  pressure  coefficients  of 
symmetric  TVD  scheme  with  Limit3  for 
projectile  with  Me,  =  1.2,  and 
Re  =  7.7xl05  using  a  70x40  C  grid. 
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Figure  7.12  The  convergence  history  of  different  TVD 
schemes  at  M.  =  1.2  using  a  70x40  grid 

(a)  upwind  TVD  scheme  with  Limit2, 

(b)  upwind  TVD  scheme  with  Limitl  and 

(c)  symmetric  TVD  scheme  with  Limit3. 
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Table  7.7  Pressure  and  Shear  Drags  for  Different  Schemes 
of  70x40  Grid  at  Ma,  =  1.2. 


Scheme 

Limiter 

Pressure 
Drag  xlO"*3 

Shear 
Drag  xlO" 

-3 

UPTVD 

Limitl 

99.51 

43.61 

UPTVD 

Limit2 

99.50 

43.52 

SYMTVD 

Limit3 

96.66 

43.51 

The  results  presented  in  the  previous  section  clearly 
illustrate  that  both  symmetric  and  upwind  TVD  schemes  are  able 
to  capture  shocks  accurately  and  efficiently.   The  numerical 
solutions  obtained  at  different  transonic  speeds  show  that  as 
the  shock  strength  increases,  the  TVD  schemes  appear  to 
produce  more  accurate  results. 


CHAPTER  VIII 
SUMMARY  AND  CONCLUDING  REMARKS 

In  the  present  dissertation,  extended  derivation  of  the 
TVD  schemes  in  three-dimensional  generalized  curvilinear 
coordinate  systems  has  been  presented.   The  derived  schemes 
are  then  simplified  and  implemented  into  a  time-dependent 
axisymmetric  thin-layer  Navier-Stokes  flow  code  by  improving 
the  treatment  of  the  numerical  dissipation  terms  to  make  the 
resulting  overall  scheme  consistent  with  the  TVD  formulas.   A 
spatially  varying  time  stepping  procedure  has  also  been 
incorporated  into  the  code  to  accelerate  the  convergence  rate 
of  the  steady  state  flow  computations.   This  improved  version 
of  the  code  was  used  to  predict  transonic  turbulent  flows  of 
Ma,  =  0.91,  0.96  and  1.2  over  an  axisymmetric  SOCBT  projectile 
at  zero  angle  of  attack.   The  numerical  results  show  good 
agreement  with  available  experimental  surface  pressure 
coefficient  data.   The  use  of  the  TVD  schemes  requires  more 
CPU  time  per  time  step  than  the  conventional  Beam-Warming 
scheme.   However,  substantial  improvements  of  the  convergence 
rate  have  been  observed  with  all  three  variants  of  TVD  scheme. 
Hence,  overall  the  TVD  schemes  are  more  economical  to  apply 
for  the  flows  studied  here.   Furthermore,  equally  accurate 
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solutions  can  be  obtained  with  the  TVD  schemes  on  less  number 
of  grid  points  compared  to  those  with  the  Beam-Warming  scheme. 
In  terms  of  the  relative  merits  among  the  three  variants  of 
the  TVD  scheme,  comparable  accuracies  are  yielded  by  both  the 
upwind  and  symmetric  schemes.   However,  it  is  found  that  for 
the  flow  with  a  weak  shock  (K^   =   0.91),  the  flux  limiter 
proposed  by  Roe  for  the  upwind  TVD  scheme  can  experience  some 
difficulties  in  convergence.  On  the  other  hand,  for  the  flow 
with  stronger  shock  (M*,  =  1.2),  the  symmetric  TVD  scheme  can 
experience  difficulties  in  the  later  stage  of  convergence.   It 
is  concluded  that  although  the  flux  limiter  proposed  by  van 
leer  for  the  upwind  TVD  scheme  may  not  be  always  the  most 
efficient  one  among  the  three  variants  studied  here,  its 
convergence  appears  less  affected  by  the  incoming  flow  Mach 
number.   The  robustness  of  this  scheme  makes  it  a  favorable 
choice  in  practical  flow  computation. 


APPENDIX  A 
NAVIER-STOKES  EQUATIONS  AND  THIN  LAYER  APPROXIMATION 

The  strong  conservation  law  form  of  the  three-dimensional 
Navier-Stokes  equations  in  Cartesian  coordinates  can  be 
written  in  the  flux  vector  form  as 


(A-l) 
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with 


125 


(A-3) 


where 


Tyy 

Tzz 
Txy 
rxz   = 

ryz  = 
£x  = 
£y   = 


=    r 


A(UX    +    Vy    +    WZ)     +    2/iUx, 

A(UX      +      Vy      +      WZ)       +      2/iVy, 

A(UX    +    Vy    +    W2)     +    2AiWz, 

yx   =   /*  (Uy   +   vx)  , 
TZX    =    M(UZ    +    Wx)  , 

rzy  =  A«(vz   +  Wy)  , 


Ur 


XX 


+    Vrxy    +    Wrxz    +    /iPr-1(7-l)"1axC2, 


UTyX     +     VTyy     +     WT  yz      +     ft  Pr"1  ( 7  - 1  )  ~1d  yC2  , 

Urzx  +  Vrzy  +  wrzz  +  ^Pr-1  (7-l)  _1a  ZC2 , 


p    is  the  density,  made  dimensionless  by  pm, 
c  is  the  local  speed  of  sound  and  c2  =  7p/p , 
u,  v,  w  are  Cartesian  velocity  components 

nondimensionalized  by  speed  of  sound  c, 
e  is  the  total  energy  nondimensionalized  by  pmc2, 
p  is  the  pressure  (  =  (7-1) [e  -  0.5p (u2+v2+w2) ] ) , 
7  is  the  ratio  of  specific  heats, 
A*  is  the  dynamics  viscosity, 
A  is  from  Stokes's  hypothesis,  =  -2/3/*, 
k    is  the  coefficient  of  thermal  conductivity, 
Re  is  the  Reynolds  number,  and 
Pr  is  the  Prandtl  number. 
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Transformed  Navier-Stokes  Equations 

To  enhance  numerical  accuracy  and  efficiency  and  to 
handle  boundary  conditions  more  easily,  the  Navier-Stokes 
equations  can  be  transformed  from  Cartesian  coordinates  to 
general  curvilinear  coordinates  where 


(A-4) 


r  =  t 

(  =   l(x,  y,  z,  t) 

v  =  r,  (x,  y,  z,  t) 

T  =  T  (x,  y,  z,  t) 


The  resulting  transformed  equations  can  be  written  in 
nondimensional  form  as 
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pV 
puV  +  ryxp 
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and 
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where  U,  V,  W  are  contravariant  velocity  components.   The 
viscous  flux  terms  are  given  by 
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where  the  components  of  the  stress  tensor  and  heat-flux  vector 
in  nondimensional  form  are  given  in  Eq. (A-3) .   The  Cartesian 
derivates  of  the  stress  tensor  can  be  expanded  in  £ ,  r? ,  and  f 
spaces  via  the  chain  rule,  for  example 


ux  =  £xu£  +  Vx^r,    +  fxuf  • 
The  metric  terms  are  defined  as 

£x  =  J(y^zr  -  yrz^)  r,x  =  J(z?yr  -  y^zr), 
£y  =  J(z,,xr  -  x^zr)  ny  =  J(x^zr  -  xzz£), 
£z  =  J(x^yr  -  y^x^)   r,z  =   J(y^xr  -  X|yr), 

(A-9)  rX    =    J(Y?Zr,     -    Z£Y„)  ?t    =    "XreX    ~    Yr^y    ?    Zr£z, 

Ty    =    J(X^Z^    -    X^Z^)  vt    =    -Xtijx    -    yrrjy    -    zrr?2, 

and  ^   =  J(^**    "  *****         rt  =  ~x'r*  "  Yrfy  -   zrrz, 

=  X£Yr/zr    +  xry?2r?    +  x^yrz^    -  x^z,,    -  x^y^zr    -  xry„Z£ 


J-l  = 


Thin-Laver  Approximation 

In  high-Reynolds-number  flows,  the  viscous  effects  are 
confined  to  a  thin  layer  near  no-slip  boundaries.   In  most 
cases,  there  are  only  enough  grid  points  to  resolve  the 
gradients  normal  to  the  body  by  clustering  the  grid  in  the 
normal  direction.   As  a  result,  the  viscous  derivatives 
associated  with  the  £  and  r,    directions  are  dropped  and  the 
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terms  in  the  f -direction  are  retained.   With  this  assumption, 
Eg. (A-5)  simplifies  to 


(A-10) 


aTq  +a^E  +a^F  +ard  =  Re_1ar§, 


where 
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where 


(A-13) 


ml  =  fx  +  fy  +  fz/ 

m2  =  fxur  +  **yvr  +  ^zwf. 

m3  =  (u2  +v2  +  w2)/2  +  Pr_1(7-1) (c2)f 


APPENDIX  B 


BEAM-WARMING  IMPLICIT  BLOCK  ADI  ALGORITHM  AND 
ARTIFICIAL  DISSIPATION 


The  finite-difference  algorithm  due  to  Beam  and  Warming 
applied  to  Eq. (A-10)  results  in  the  following  approximate 
factorization : 

(B-l)      Le  L,,  Lr  Aq  =  (RHS)n, 
where 

L^  =  (I  +  hs^kn   +  DjJ^), 


L„  =  (I  +  hS„Bn  +  Di|„)f 

Lr  =  (I  +  h*rCn  -  hRe-1*rJ-1MnJ  +  Dilr), 
(RHS)n  =  -At(^En  +  ^Fn  +  5rdn  -Re-ifij-finj  -  De, 


and 


where  6    is  the  central-difference  operator  and  A  and  v  are 
forward  and  backward-difference  operators,  and  h  =  At 
corresponds  to  the  first-order  time-accurate  Euler  implicit 
scheme.   DjJ^,  Di  |  n    and  Di  |  j-  are  the  second-order  implicit  and 
De  is  the  fourth-order  explicit  smoothing  operator. 

The  Jacobian  matrices  An,  Bn,  Cn  and  Mn  are  obtained  by 
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linearizing  the  flux  vectors  £n,  $n   &n  and  iln  in  time  such 
that 
(B-2) 

fcn+1  =  £n  +  £nA£  +  0(h2), 

f.n+1  .  fn  +  fcnA£  +  0(h2) , 

£n+l  =  £n  +  £nA*  +  0(h2), 

Re-l§n+l  =  Re-l[gn  +  j-lftn^nj  +  0(h2)  , 

where  k  =  dt/dq,    6  =  d$/dq,    £  =  dt/dq,    and  ft  =  3§/dq  are  the 
flux  Jacobians. 

The  Jacobian  matrices  are  i  or  I  or  6  = 
(B-3) 
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with  k   =  £,  or  i)  or  f  for  A,  §,  or,  6  respectively, 

The  viscous  flux  Jacobian  is 
(B-4) 
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where 


m21  =  -ttlar  (u/p)-a2af  (v/p)-a3df  (w/p)  , 
m31    =    -<*2af  (u/p)-a43r  (v/p)-a53r  (w/p)  , 
m41    =    -«3*J (u/p)-a53r (v/p)-a6af (w/p ) , 
m44    =   a43f  (p-1)  , 
m51  =  "O^f [-(e/p2)+(u2+v2+w2)/p] 

-aiar (u2/p)-a4ar (v2/p)-a6ar (w2/p) 

-2a2af (uv/p)-2a3ar (uw/p)-2a5ar (vw/p ) , 
m52    ■    -a0af (u/p)     "   ^21/ 

m53  =  -«oa$-(v/p)    -  m31, 
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<*3  =  (^/3)rxfz# 

"  y 


•4    =   ^tTx2    +     (V3)fy2    +    fz2], 


<»5  =  (M/3)rxrzf 

«6   =  M[fx2    ♦    fy2    +    (V3)rz2]. 

Artificial  Dissipation 

The  numerical  solution  of  the  Euler/Navier-Stokes 
equations  for  flows  involving  captured  shock  waves  require  the 
use  of  numerical  or  artificial  dissipation.   The  reason  for 
adding  artifical  dissipation  is  to  control  strong  nonlinear 
effects  such  as  shocks. 

The  block  algorithm  which  was  given  in  Eq. (B-l)  uses  a 
constant  coefficient  dissipation  model.   In  this  approach,  an 
explicit  fourth-order  dissipation 

(B-5)  De  qn  =  £e  At  ^[(V^)2  +  (V^A,,)2  +  (VrAr)2]  Jqn, 

is  added  to  the  right-hand  side  of  the  equation  and  an 
implicit  second-order  dissipation 

(B-6)       Di|k  =  -eij-l  AkvkJ,     *  -  f ,*,t. 

is  inserted  into  the  respective  implicit  block  operators. 
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